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Abstract 

A general-purpose, self-adapting Monte Carlo (MC) event generator (simulator) 
Foam is described. The high efficiency of the MC, that is small maximum weight or 
variance of the MC weight is achieved by means of dividing the integration domain 
into small cells. The cells can be n-dimensional simplices, hyperrectangles or a 
Cartesian product of them. The grid of cells, called "foam", is produced in the 
process of the binary split of the cells. The choice of the next cell to be divided and 
the position/direction of the division hyperplane is driven by the algorithm which 
optimizes the ratio of the maximum weight to the average weight or (optionally) the 
total variance. The algorithm is able to deal, in principle, with an arbitrary pattern 
of the singularities in the distribution. As any MC generator. Foam can also be used 
for the MC integration. With the typical personal computer CPU, the program 
is able to perform adaptive integration/simulation of a relatively small number 
of dimensions (< 16). With the continuing progress in CPU power, this limit will 
inevitably get shifted to ever higher dimensions. Foam program is aimed (and already 
tested) as a component of the MC event generators for the high energy physics 
experiments. A few simple examples of the related applications are presented. Foam 
code is written in fully object-oriented style, in the C++ language. Two other 
versions with a slightly limited functionality, are available in the Fortran?? language. 
The source codes are available from http://jadach.home.cern.ch/jadach/. 
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PROGRAM SUMMARY 



Title of the program: 
Foam, version 2.05. 
Computer: 

any computer with C++ or Fortran 77 compilers and UNIX operating system 
Operating system: 

UNIX; program was tested under Linux 6.x. 
Programming languages used: 

ANSI C++ and FORTRAN 77 with popular extensions such as long names, long code 
lines, etc. 

High-speed storage required: 
< 50 MB 

No. of hnes in combined program and test deck: 
4235 lines of C++ code and 9826 hnes of F77 code. 
Keywords: 

Monte Carlo (MC) simulation and generation, particle physics, phase space. 

Nature of the physical problem: 

Monte Carlo simulation or generation of unweighted (weight equal 1) events is a standard 
problem in many areas of research. It is highly desirable to have in the program library 
a general-purpose numerical tool (program) with a MC generation algorithm featuring 
built-in capability of adjusting automatically the generation procedure to an arbitrary 
pattern of singularities in the probability distribution. Our primary goal is simulation 
of the differential distribution in the multiparticle Lorentz invariant phase space for the 
purpose of comparison between Quantum Field Theory prediction and experiments in the 
high energy experiments. 
Method of solution: 

In the algorithm, a grid (foam) of cells is built in the process of the binary split of the 
cells. The resulting foam is adapted automatically to the shape of the integrand in such a 
way that the resulting ratio of average weight to maximum weight or variance to average 
weight is arbitrarily good. The above algorithm, is a substantial improvement on the 
previous version in Ref . [1] . The division of the cell is improved and, in addition to cells 
of a simplical shape of Ref. [1] , a hyperrectangular cells are also used. 
Restrictions on the complexity of the problem: 

The program is memory- hungry and therefore limited, at present, to relatively small 
dimensions < 16. (In Foam 1.x of Ref. [1] the dimension was limited to n < 6.) 
Typical running time: 

The CPU time necessary to build up a foam of cells depends strongly on the number of 
dimensions and the requested size of the grid. On the PC with a 550 MHz Intel chip, 
it takes about 30 seconds to build a hyperrectangular grid of 10000 cells for a simple 
3-dimensional distribution. 

[1] S. Jadach, Comput. Phys. Commun. 130, 244 (2000). 
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1 Introduction 



This work describe a new version of an algorithm for producing random points according 
to an arbitrary, user-defined distribution in the n-dimensional space - much improved 
with respect to the original version of Ref. The new implementation is realized in the 
C++ programming language in a fully object-oriented mannerQ. Since the changes, with 
respect to Ref. [Q, both in the algorithm and in the implementation are quite essential, 
a complete description of the method and the new code is provided, instead of only an 
update of Ref. IQ . 

For the problem of function minimization, integration (quadrature), there are plenty 
of general-purpose programs, which can be applied to an arbitrary user-defined function. 
"General-purpose" means that all these tools work, in principle, for a very wide range of 
user-defined functions. For the multidimensional Monte Carlo simulation problem, that 
is for the problem of randomly generating points according to a given n-dimensional dis- 
tribution, there is precious few examples of the General- Purpose Monte Carlo Simulators 
(GPMCS), that is programs that work (in principle) for arbitrary distribution As 
in these works, here we are concerned mainly with the MC applications to high energy 
physics, that is to simulation of the differential distributions in the multi-particle Lorentz 
invariant phase space provided the Quantum Field Theory, see also classic Ref. ^ on this 
subject. An example of the work on GPMCS applied to other fields see the interesting 
worksQ of Refs. g^. 



Let us also note that the two-dimensional cellular MC sampler VESK02 with the 



primitive binary split was already included in the program LESKO-F of Ref. W^ long 



time ago. Very similar programs were also described in Refs. ||12| , |T3|0. Still another 



very interesing class of GPMCSs for the high energy physics, based of the Metropolis 



algorith |14], is described in Ref. |15]. 



Two essential reasons for a realtive scarcity of GPMCSs, as compared to other pro- 
gramming tools, are the lack of novel ideas for an efficient algorithm and the need of much 
CPU power and memory - only recently available or affordable. 

GPMCS is essentially a random-number generator of points in multidimensional space 
with a non-uniform user-defined probability distribution. Inevitably the GPMCS works 
in two stages: exploration and generatior^. During an exploration phase, the GPMCS 
is "digesting" the entire shape of the ra-dimensional distribution p{x) to be generated, 
memorizing it as efficiently as possible, using all available CPU processing power and 
memory^ In Foam, the exploration phase is the phase of the build-up of the system of 

^ The early C++ version of Fosun was coded by M. Ciesla and M. Slusarczyk It was a translation 
of a version 1.x from Fortran?? to C++. Analogous translation to the JAVA language was also done. 

^In these works there is far more emphasis on the parallel computing aspects of the integration (not 
simulation) than on the cell geometry, as compared with our work. 

•^We thank S. Kawabata for drawing our attention to these works. 

"^Exploration and generation could be done simultaneously, at the expense of complications in the 
algorithm and the code. 

^The procedure of memorizing multidimensional distribution p{x) > is a kind of interpolation, in 
which the grid of cells is denser in places where the distribution peaks and/or varies strongly. 
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cells covering entirely the integration space, which will be called "foam", produced in 
the process of the binary split of the cells. In the generation phase, GPMCS provides 
a method of the MC generation of the random points x exactly according to p{x). The 
vector X = X = {xi,X2, ■ ■ . , Xn) will also be called in the following a Monte Carlo event. 
In Foam, the MC generation is very simple: a cell is chosen randomly, and next, a point 
is generated within the cell with uniform distribution; see below for more details. The 
value of the integrand is already estimated in the exploration; it can be calculated with 
an arbitrary precision in the generation phase. 

During the exploration Foam constructs a distribution p'{x), which is uniform within 
each cell, and is used for the MC generation. Events are weighted with the weight w = 
p/p'. The quality of the distribution of this weight, measured in terms of the weight 
distribution parameters, such as variance and ratio of maximum to average, is determined 
by the quality of the exploration. The basic principle of the Foam algorithm is that the 
parameters of the anticipated "target weight distribution" in the MC generation phase are 
used as a driving force guiding the cell build-up (exploration). In the case of a successful 
exploration, weighted MC events can be turned efficiently into unweighted ones with the 
usual rejection method, that is with a small rejection rate. 

Since the exploration phase may be CPU-time consuming, it is natural to expect that 
GPMCS has a built-in mechanism of persistency, i.e., there is a mechanism of writing 
into a mass storage (computer disk) the whole information on the memorized shape of 
the distribution obtained from the exploration phase, such that the generation of the 
MC events can be (re)started at any later time, without any need of repeating the time- 
consuming exploration. One small step further is to require that the generation of events 
with GPMCS can be stopped at any time, the entire status of the GPMCS can be written 
on the disk, and the generation of the next event can be resumed at any later time upon 
reading the stored information; the next generated event will be such as if there had been 
no break in the generation process. In fact, this is what we shall really mean in the 
following as a persistency mechanism for GPMCS, and what is actually implemented in 
the Foam. There, the persistency is realized using the ROOTQ package p6[| . 

The GPMCS programs will always be limited to "small dimensions" . With currently 
available computers, "small" means in practice n < 10, up to n < 16 for "mildly" singular 
distributions. This is already quite satisfactory, especially if we remember that this limit 
will pushed higher, as the available hardware gets more powerful, without any need of 
modifying the existing code|]. Twelve years from now, with portable computer featuring 
a 100 GHz processor and 1 TByte disk the same version of Foam will work efficiently for 
even higher dimensions. 

Foam has been developed having in mind that it will be used as part of a bigger MC 
program, typically, to generate a subset of variables in which a model distribution is the 
most singular (has strong peaks). This is why we are not so much concerned by the 

^ The use of ROOT is optional in Foam. However, version of Foam without ROOT docs not feature 
any kind of persistency. 

^The present implementation of Foam is fully based on the dynamic allocation of the memory and the 
space dimension is a user-defined parameter. 
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fact that the cellular algorithm of Foam is inefficient for, say, 150 variables. The user is 
supposed to select n < 16 "wild variables" and apply Foam to them. For the remaining 
"mild variables" , Foam may merely serve as a provider of the uniform random numbers, 
if the user of Foam wishes to exploit that option. On the other hand, for smaller MC 
problems. Foam may play the role of a "stand-alone MC generator" or "stand-alone MC 
integrator". Also, from the following description of the various modes of the use of Foam 
it will be clear that the subprogram providing the model distribution to Foam can have 
quite a complicated structure. Nevertheless, this user-provided part of the program will be 
smaller than a solution without Foam, because Foam provides for essential functionalities 
concerning the optimization of the MC weight distribution. This remark is especially true 
for the case of implementation of multibranching with the help of Foam. 

It is worth mentioning that Foam is not based on the "principle of factorizability" of the 
integrand distribution, p(x) = Hi Pj(^i); which VEGAS-family programs are built 

The outline of the paper is the following: in Section 2 we describe the cellular Foam 
algorithm, delegating the description of the cell division procedure to Section 3. Section 4 
is devoted to the description of the Foam code in C-I--I-. The use of Foam is described in 
Section 5 and examples of the numerical results (MC efficiency) are given in Section 6. 
Conclusions and Appendices on the variance minimization finalize the paper. 

2 The Foam algorithm 

As already mentioned, the execution of the Foam algorithm is clearly separated into the 
first stage of the "distribution exploration" consisting of the "build-up of the foam of 
cells", which in a sense memorizes the n- dimensional shape of the distribution, and the 
second stage of the actual "MC generation", see Fig. 0. The most essential part of the 
present Foam algorithm is the procedure of the binary split of the cell, in which it is 
decided which cell is picked up for the next split and the necessary parameters of the 
geometry of the cell split are determined. This part of the Foam algorithm description 
is delegated to the next section. In the present section we describe other, more general, 
aspects of the Foam algorithm. 

2.1 Cellular exploration of the distribution 

The most obvious method to minimize the variance or maximum weight of the target 
weight distribution in the MC generation, already considered some 40 years ago, is to split 
the integration domain into many cells, so that the distribution p{x) is approximated by 
p'{x), which is constant within each cellQ. This is a cellular class of the general-purpose 
MC algorithms^. 

® For the MC simulation, our main aim, a more sophisticated interpolation of p{x) within a cell does 
not seem to be worth the effort - it would be interesting if our main aim was the integration of p{x). 

^The term "stratified sampling" , used in the literature, has in our opinion a narrower meaning than 
"cellular class" . 
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I Build-up of the foam of cells 

Split root cell if necessary 

Choose next cell for the split 

i 

MC exploration of the cell 

Generate series of MC events inside a cell 
Choose best direction (division edge) 
Find out best division ratio (division plane) 



i- , 1 

Generate MC event 

Choose randomly a cell | 
Choose randomly a point inside a cell ! 

I 'J 

Figure 1: Two stages in the cellular algorithm of Foam. 

The immediate questions are: What kind or shapes of the cells to use and how to 
cover the integration domain with cells? The reader may find in Ref. an example of a 
rather general discussion of these questions. In the Foam program the user may opt for one 
of the three geometries of the cells: (1) simplices, (2) hyperrectangles and (3) Cartesian 
products of simplices and hyperrectangles. For these particular types of cells there exists 
an efficient method of parametrizing them in the computer memory and handling their 
geometry. 

The system of many cells can be created and reorganized all at once, as in VEGAS-type 
programs or in a more evolutionary way, as the cell split process of this work. In the 

Foam algorithm we rely on the binary split of cells. Starting from the entire integration 
domain (unit hyperrectangle or simplex) cells are split into two daughter cells, step by 
step, until the user-defined memory limit is reached. The choice of next cell to be split 
and the geometry of the split in the exploration phase are driven by the "target weight 
distribution" of the generation process, see Section ^. The important advantage of any cell 
split algorithm is that it assures automatically the full coverage of the integration domain 
- simply because the primary root cell is identical to the entire integration domain and 
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the two daughter cells always cover the parent cell entirely. The problem of blind spots 
discussed in Ref. is avoided by construction. 

In the early version of Foam there was a possibility in the algorithm that the 
"unsuccessful" branch in the tree of all cells can be erased and rebuild. This was called 
"collapse" and "rebuild" . In the present version this option was removecf^ because the 
experience with many testing functions has shown that the algorithm of the cell build-up 
is rather "deterministic" and the "rebuild" procedure was usually leading to a new branch 
of foam with about the same features as the old one. 

Let us finally remark that the version of the cellular algorithm presented in this paper 
is, in fact, result of many experiments with different variants of the cellular algorithm. 
The presented version is the best one out of several development versions. In the code 
one may still see some "hooks" and unused features (class members or methods) related 
to these alternative variants. We have left them just in the case that some new idea of 
improving the algorithm emerges, or for certain kinds of debugging/testing. 

2.2 Variance reduction versus maximum weight reduction 

In the construction of the Foam algorithm, most effort was invested into a minimization 
of the ratio of the maximum weight to the average weight Wmax/ (w). This parameter is 
essential, if we want to transform variable-weight events into w = 1 events, at the latter 
stage of the MC generation|^. 

Minimizing the maximum weight Wmax is not the same as minimizing the variance a = 
a/ {w'^) — (w)'^. Minimizing Wmax may be more difficult - but it is worth an effort because 
Foam is really meant to be a part of a bigger MC program, where it is usually essential that 
the "inner part" of the program provides events with an excellent weight distribution, or 
even w = 1 events. Nevertheless, minimizing the variance is also implemented in Foam 
and available optionally. It can be useful if one is satisfied with the variable-weight events, 
and/or if the main aim is the evaluation of the integral and not the MC simulation. 

The difference between the above two options is well illustrated in Fig. ^ which shows 
two examples of the evolution of the MC weight distribution due to a gradual increase 
of the number of cells. For the default configuration. Foam optimises the ratio Wmax/ (w) . 
This case is shown in plots (a-c) in Fig. |^. Here, the weight distribution features a sharper 
and sharper drop of the weight distribution at u; = 1, while increasing the number of cells. 
Also, the average weight increases gradually and the weight distribution gets narrower. 
The optional case of the optimization of (t/{w) is shown in plots (d-f) of Fig. |^. In this 
case the variance is decreasing with the growing numbers of cells. On the other hand, the 
maximum weight is noticably higher than before. All weight distributions were obtained 
for the same 2-dimensional testing function Pb{x), used also in Section |^. 



^"^A "flush method" , which erases the entire foam of cells from the computer memory and allows for 
its reinitialization is, however, available. 

^^We provide optionally in Foam for the rejection leading to w = 1 events. 
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Figure 2: Weight distribution of the Foam for the default option with the maximum weight optimiza- 
tion (a-c) compared to analogous distributions obtained for an option with the variance optimization 
(d-f). Number of cells is 200, 2000 and 20000 for (a-c) and (d-f), correspondingly. 



2.3 Hyperrectangles or simplices? 

In Ref. simplical cells have been chosen instead of simpler hyperrectangles, mainly 
because of the author's "prejudice" that simplices may adapt more efficiently to com- 
phcated singularities in the distribution p(x) spanned along subspaces, not necessarily 
parallel to the axes of the global reference frame. Hyperrectangles tend to remember the 
orientation of the parent hyperrectangle, while simplices feature, in principle, a kind of 
"angular mobility", i.e. they can forget the orientation of great-grand-parents, and adapt 
to the orientation of the singularity in p{x). However, an experience with tens of testing 
functions has shown that in many cases hyperrectangles provide the same or even better 
final MC efficiency than simplices, for the same number of cells. Moreover, simplices have 
certain additional disadvantages. At present. Foam with simplices is practically limited 
to rather low dimensions n < 5, because in most cases the entire integration domain is a 
unit hyperrectangle, which has to be divided into n\ simplices, where n! quickly becomes 
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a large number0. This limitation is, of course, not valid, if the entire integration domain 
is actually a simplex of high dimensionality instead of a hyperrectangle. (Foam can be 
configured to start cell evolution from a simplex or Cartesian product of a simplex and 
a hyperrectangle.) Furthermore, geometry manipulations in the simplical case require 
calculation of many determinants - this slows down the program execution at higher 
dimensions. In addition, in the present implementation, the memory consumption in a 
simplical foam build-up is ~ 16 x n Bytes/cell, while for hyperrectangles we have found a 
method that limits memory consumption to below ~ 80 Bytes/cell independently of n, see 
Section p.6| . We can therefore reach easily the level of 10^ hyperrectangular cells at any 
dimension (in practice n < 16) and about 50, 000 simplical cells, for n < 5. As we see, 
hyperrectangular foam seems to win on many fronts. Nevertheless, we keep simplical foam 
as an option, because in certain applications one will definitely encounter distributions 
for which it turns out to be more efficient to use simplices, in spite of all their limitations, 
at least for a subset of the integration variables. 



2.4 Build up of the foam and data organization 

The foam of cells is built-up by starting from the root cell, which is the entire integration 
domain, through a process of binary split of a parent cell into two daughter cells. The 
root cell is either a unit hypercube < < 1 (default) or a simplex < Xi < X2 < 2:3 < 
• • • < < 1. Also a Cartesian product of these two shapes available on option. Any 
cell being a product of the cell split can be also a hyperrectangle, a simplex or Cartesian 
product of the fc-dimensional hyperrectangle and n-dimensional simplex, with the total 
dimensionality k + n. If the starting root cell is a hypercube and cells are simplical (or 
of mixed type) then the root cell is immediately divided into n! simplical (or mixed type) 
cells. 

Each cell is explored immediately after its creation. In the exploration of the cell, 
about 100-1000 MC events (the user may define this number) are generated inside the cell 
with flat (uniform) distribution and using MC weight equal to p{x); certain averages and 
certain integrals over the cell are estimated. Also, the best geometry of the binary split of 
the cell is established and recorded for future use. In this way, every created cell is ready 
for an immediate split. The determination of the best split is described in fine detail in 
Section ^ In the exploration the estimate of the integral Rj = f^^ p{x)dx"' is calculated 
for each cell uj. Far more important is another functional Rioss\i = f^^ Pioss{x)dx"', see 
Section ^ for its definition, which determines the evolution of the foam and the split of the 
cell. Next cell to be divided into two is a cell chosen randomly, according to a probability 
proportional to Rioss\i or, optionally, a cell with the largest Rioss\i- 

The process of cell division continues until the user-defined maximum number A^^ of 
cells is reached. A^ includes, in addition to the normal cells called active cell, also all cells 
that were split, i.e. all parents, grand-parents and great-grand-parents inluding root cell. 

Mapping of the hyperrectangle into a simplex is possible, however, one should use transformation 
with \dx/dy\ ^ 0, in order to avoid producing nasty singularities located at the vertices, edges and walls 
of the simplex. 
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which we shall call inactive cells. Usually, when we refer to cells, we mean both active and 
inactive ones. Keeping inactive cells in the record may look like a waste of the memory, 
but owing to the binary character of the cell split, the loss is only a mere factor of 2. 
There are many reasons, which make it profitable to keep all inactive cells (including the 
root cell) in the memory; in particular, as we shall see in Section |2.6| , keeping all cells 
in the record will help us encode all cells in the memory in an economic way, that at 
higher dimensions we finally gain in terms of total consumption of the CPU memory. 
Furthermore, for certain quantities that are the integrals over the cell, such as Rj, we 
do the following: just after the split, when a new, more precise value of Rj is known for 
the daughter cells - the value of the Rj of the parent cell is updated with the sum of 
the contributions from two daughter cells. This correcting procedure is repeated for all 
ancestor cells up to the root cell. In this way, the root cell (and any other inactive cell) 
always keeps track of the actual value of the total Rj during the whole foam build-up 
process. This can be done for any other integral quantity as well, and can be exploited 
for various purposes. 

Since the maximum number of cells Nc is defined at the beginning of the foam build-up, 
all the cell objects and/or other related objects (vertices) are allocated in the computer 
memory at once, at the very beginning of the cell build-up. On the other hand, the 
cell objects are organized as a multiply-linked list, with pointers towards parents and 
daughters. In addition, an array of pointers to all active cells is created at the end of the 
foam build-up, in order to speed-up the MC generation. 

Let us now briefly explain how the geometry of an individual cell is parametrized and 
stored in the memory. It is relatively easy to parametrize an n-dimensional hyperrectangle 
or simplex in a way that does not require much computer memory. An n-dimensional 
simplex is fully determined by its n + 1 vertices. Since most of the vertices are shared 
by the adjacent simplices, the most efficient method is to build an array of all vertices^ 
Vk,K = 1,2, Ny, each of them being an ?7,-component vector, and to define every 
simplex as a list of n + 1 vertex indices (integers or pointers) Ki,K2, ...,Kn+i- For N^. 
simplical cells resulting from the binary split of a single "root" simplex cell, the number 
of vertices is n + 1 + N^, because each binary split adds one new vertex. (We include in 
Nc also cells that have got split). The interior points of the simplex are parametrized as 
follows 

n 

S=Y,HVk,-Vk,), A,>0, J]A,<1, 1 = 1,2,. ..,n, (1) 

using basis vectors relative to the p-th vertex. The above method would be inefficient 
for n-dimensional hyperrectangles, because memorizing all 2" vertices would require too 
much memory at higher dimensions. Instead, we use another way of parametrization: each 
hyperrectangle is defined by the n-dimensional vector q defining the origin of the cell and 
another vector h = [hi, h2, hn), where each component is the length of the hyper- 
rectangle along the i-th direction. This is even clearer from the explicit parametrization 

^■^In the foam build-up every new vertex is appended at the end of the array. 
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of the interior of the hyperrectangle: 

Xi = qi + \ihi, < Ai < 1, i = l,2,...,k. (2) 

For cells with mixed topology, we apply Eq. (0) for i = 1,2, ... ,k and Eq. (P for 
i = k + 1, k + 2, . . . , k + n. In Section we describe an optional method of storing 
hyperrectangular cells, in which just two integer numbers are recorded instead of two vec- 
tors q and h (two 2-Byte integers instead of 2n of 8-Byte floating-point numbers). This 
method is implemented for the hyperrectangular part of the space only. 

2.5 Monte Carlo generation 

Once the build-up of the cells is finished, the Monte Carlo generation takes place. There 
is no need for any reorganization of the cells. MC generation can be started immediately. 
The only thing, that is done at the very end of the foam build-up is the preparation of 
the list of pointers to all active cells and the array of corresponding R'j. 

The MC point is generated in two steps. First, a cell is chosen with a probability 
proportional to R'j = J^^q^Hj p'{x) and next a MC point x is chosen with uniform prob- 
ability inside the cell. The MC weight w = p'{x)/p{x) is associated with the event. For 
a successful foam of cells, the MC weight is close to 1 and the user may turn weighted 
events into w = 1 event by means of the rejection method (with the acceptance rate 
~ {w) / Wma.x) ■ Foam can do this also for the user. However, the user can sometimes orga- 
nize the calculation of the (w) and bookkeeping of other parameters of the weight better, 
in a way that best fits his own aims. This is why the mode of variable weights MC events 
is also available. The total integral, usually necessary for the proper normalization of the 
MC sample, is calculated using R = R'{w). Foam program provides both the exact value 
of the R' and the MC estimate of the integral R. 

2.6 Economic use of the computer memory 

The actual implementation of the single cell object occupies about 80 Bytes (it could be 
shrinked to about 40 Bytes, if necessary) of various integer and double precision attributes, 
plus the dimension-dependent part. In the case of a simplical cell, each new cell adds one 
n-component double-precision vector (vertex) and the total memory consumption is there- 
fore of order (80 + 8 x n) Bytes/cell. For n = 5 and lOOK cells it is therefore ~ 15 MBytes 
of the memory, still an affordable amount. For the hyperrectangle cells we have to count 
two ra-component double-precision vectors per cell, that is about (80 + 16 x n) Bytes/cell. 
For the 10^ cells and n = 15 that would mean ~ 340 MBytes for the entire foam of cells, 
which could be annoying. Fortunately, we have found a method of substantially reducing 
the memory consumption for a hyperrectangular foam. As discussed in Section |] the 
geometry of the cell division is fully determined in terms of two integers: one of them 
is the index of an edge to which the division plane is perpendicular and the other one 
defines the position of a division plane. The position parameter is a rational number, and 
only the integer numerator has to be remembered, while the denominator is common to 
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all cells. The above two integers define uniquely the position of the two daughter cells 
relative to a parent cell. With this method the memory consumption is down to about 
80 Bytes/cell, independently of n in the present implement at iorf^ There is, however, a 
price to be payed in terms of CPU time. For the generation of the point inside a cell, or 
even the evaluation of the weight, we need the "absolute" components of x, that is in the 
reference frame of the root cell, not relative to vertices of the cell. It is therefore necessary 
to use a procedure (a method in the class of cells) to construct the absolute position of a 
given cell "in flight" . This is done by means of tracing all ancestors of a given cell up to 
the root cell and translating position and size with respect to its parent into absolute ones, 
step by step, finally relative to the root cell. It is implemented exploiting the organization 
of cell objects into a linked binary tree. The average number of ancestor cells to be traced 
back from a given active cell up to the root cell for Nc = 10^ cells is on the average about 
ln2 Nc ~ 20. This may cause about 20% increase in the CPU time of the MC generation 
- an affordable price, if we remember that the MC efficiency improves mainly with the 
increase of Nc. In principle, this kind of memory-saving arrangement is also possible for 
simplical cells; however, in this case the CPU time overhead would be larger, because of 
the necessity of the full linear transformation at each step, on the way from a given cell up 
to the root cell. In the case of hyperrectangular cells the transformation is much simpler 
(and faster); it is the translation and/or dilatation along a single spatial direction at each 
step. 




Figure 3: Inhibited cell division for first variable, that is for xi (right). Foam with 250 cells. 



*In fact, it can be reduced below 40 Bytes/cell, if really necessary. 
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2.7 CPU time saving solution 



The final MC efficiency is improved mainly by increasing the number of cells Nc- The 
CPU time of the cell build-up is T ~ n x A^^^ x Ngamp, where N^amp is the number of 
MC events used in the exploration of each newly created cell. The important practical 
question is: Can one somehow reduce Nsamp without much loss of final MC efficiency, in 
order to be able to increase Nc, within the same CPU time budget? 

A simple solution is the following: during the MC exploration of a new cell we con- 
tinuously monitor an accumulated "number of effective events with w = 1" defined 
as Neff = i^Wi)^/ '^wf, and terminate cell exploration whenQ Nf^jf/riun > 25, where 
nbin is the number of bins in each histogram, which is used to estimate the best division 
direction/edge parameters. This method helps to cut on the total CPU time, because the 
increase of Nsamp is not wasted for cells in which the distribution p[x) is already varying 
very little. At the later stage of the foam evolution this happens quite often. In this 
method the user may set Ngamp to a very high value and the program will distribute the 
total CPU time (in terms of Ngamp) among all cells economically, giving more CPU time 
to those cells really need it, i.e. to cells with the stronger variation of p{x). 

2.8 Inhibited variables — flat dependence 

In some cases the user may not want Foam to intervene into certain variables in the 
distribution p{x), simply because there is little or no dependence on them in p{x). The 
user may draw, of course, these variables directly from any uniform random number 
generator. He may, however, find it more convenient to get them from the Foam program. 
This is easily implemented in Foam: any variable Xi may be "inhibited" for the purpose 
of cell-splitting procedure. In the Foam code it is actually done in such a way that Foam 
excludes this variable (edge) from the procedure of determining the best binary division 
of the cell. This provision makes practical sense mainly for the hyperrectangular part of 
the variable subspace. 

In Fig. ^ we show two 2-dimensional foams (250 cells) for the same testing distribution 
p{x) (two Gaussian peaks on the diagonal). In one of them (right plot) we have inhibited 
a split in the ffist variable, that is for xi. 

2.9 Predefined split points — provision for very narrow peaks 

In the practical applications, see Refs. [|T7|,p^, one may encounter in certain variables 
extremely narrow spikes (narrow resonances). The Foam exploration algorithm may find 
it difficult to locate these spikes with the usual method of MC sampling in the cells, at the 
early stage of the foam build-up. For very narrow spikes, or too low number of requested 
cells, it may not find them at all! The user usually knows in advance the position of these 
spikes and the Foam should have a built-in mechanism to exploit this knowledge. 

^^The actual limit of equivalent events per bin is the user-defined parameter, not necessarily equal 
to 25. 
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Figure 4: Predefined division points at xi = 0.30,0.40, and 0.65, for 2000 cells. 



The solution is simple. (It applies for the hyperrectangular subspace of the parameter 
space only.) The user has a possibility to provide Foam, for each variable, with the list of 
a number of predefined values: the first splitting positions of the root cell. In the Foam 
algorithm, it is checked if the list of predefined division points is not empty. If it is the 
case, then instead of adopting the division parameter from the usual procedure described 
in Section ^, Foam takes the division parameter from the list, and removes it from the 
list. In this way the first few division points are taken from the "user-defined menu", 
if available, and the next ones are chosen with the usual methods. For narrow spikes 
this method helps Foam to locate and surround them with as a dense group of cells as 
necessary. 

In Fig. ^ we show an example with two Gaussian peaks in which we requested the 
Foam program to use the three predefined division points for the Xi variable. They are 
clearly seen as three vertical division lines dividing the entire root cell. In the present 
case, peaks are not so narrow and there is no real need for a predefined division. The 
example is just to illustrate the principle of the method. 
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2.10 Mapping of variables 

If the structure of the singularities is known and/or Foam is unable to get a reasonable 
weight distribution for a reasonable number of cells, it is then worth performing an ad- 
ditional change of variables, so that the transformation Jacobian compensates for the 
singularities, at least partly. In such a case the user subprogram provides Foam with the 
distribution 



dx'^^'> (y) 



dy 



(3) 



instead of the original p{x) = d'"' p / dx"" . For each vector y generated by Foam, the image 
vector X is well known in the subprogram calculating p*{y). A mechanism for exporting 
X to the outside world usually has to be provided by the user, because Foam cannot help 
- it does not know anything about x; it only knows y. 

Note that in the limiting case of the "ideal mapping" we have 



dx{y) 



dy 



" (4) 



p{x 



consequently, p*{y) = R and in this case Foam would play merely the role of a provider of 
the random numbers for y. 

The user of Foam may also need to apply mapping in the case of a "weak" integrable 
singularity in the distribution p like log(x) or ^/x. Foam can deal with them by brute 
force, at the expense of a larger number of cells. However, a wiser approach is to apply 
mapping, in order to remove such singularities from the distribution. 

In the next section we shall describe how to combine the mapping method with the 
multibranching. Such a mixture is well known as the most powerful method of improving 
the efficiency of the Monte Carlo method. 

2.11 Provisions for the multibranching 

In the following we elaborate on the various methods of implementing multibranching MC 
method p!9|-p2| with the help of Foam. 



2.11.1 Single discrete variable 

As a warm-up exercise, let us consider the question: Is Foam capable to generate (and 
sum up) a discrete variable i = 1,2, . . . , N, according to the (unnormalized) distribution 
Ti, . . . ,ri\f7 Of course it can. The simplest way is to define an auxiliary 1-dimensional 
distribution 

p{x) = ri, for ^ <x< ^, z = l,2,...,Ar. (5) 

The user subprogram providing the above p{x) is trivial. If plotted, this p{x) would 
look like a histogram with equal-width bins. Foam will build up its own grid of cells 



16 



(intervals), and if we request a high enough number of cells (that is Nc > N), it will 
approximate the above p{x) very well, with its own "histogram-like" distribution p'{x). 
However, the Foam approximation will never be ideal, because Foam is not able to detect 
the exact position of the discontinuities in p{x). (Nevertheless, this will be a workable 
solution with a very good weight distribution.) The present Foam algorithm provides for an 
essential improvement: one may predefine the division points as x^*-* = i/N, i = 1, . . . , N, 
and set the number of cells to be Nc > N. In such a case Foam will define its cells matching 
exactly the shape of p{x). It will generate points with w = 1 and provide the exact sum 
R = R' = ^Ti, already at the end of the foam build-up. In the MC generation. Foam 
will randomly generate variable x G (0, 1), which can be translated easily into discrete i 
l<i< N. 



2.11.2 Discrete and continuous variables 

How does the above extend to the case of the distribution p depending on one discrete 
variable and the usual n continuous variables? For such a distribution p{yi, . . . ,yn,i) we 
define 

t — 1 t 

p{xi,X2,...,Xn+i) = p{xi,...,Xn,i), for < a;„,+i < — , z = 1, 2, . . . , A^, (6) 

in completely analogy with Eq. (^. As previously, we provide for the variable Xn+i 
a list of predefined division points = i/N, i = 1,2, ...,A^, and, of course, we 

request Nc » N. There is still one small problem: Foam may "by mistake" perform 
an unnecessary cell division for the variable Xn+i, simply due to statistical errors in the 
"projection histogram" described in Section p.4.1| . This problem is solved in Foam in 



an elegant way: in addition to providing for Xn+i predefined division points, the user 



of Foam may declare Xn+i as an "inhibited variable" in the sense of Section |2]^. In 
this case Foam will still split cells according to a list of predefined division points for 
Xn+i, but will not perform any additional division in this variable! (The translation of the 
continuous Xn+i to the discrete i is done as usual.) The above method is the basic method 
of implementation of the "multibranching" (or "multichannel" ) MC method using Foam. 
Let us call it "predefined and inhibited division" , for short PAID method. We shall also 
describe below how to combine the PAID method with mapping, etc. 

In order to appreciate more fully the advantages of PAID, let us consider a more 
straightforward implementation of the multibranching. In the object-oriented environ- 
ment one may easily construct instances of the Foam object, each of them for the 
n-dimensional function p{xi, . . . , x„, i), initialize them (creating a separate foam of cells) 
and generate event {x,i) with the associated weight Wi{x). Index i can be chosen accord- 
ing to probability pi = R'J R'j, where R[ are provided by the i-th object of the Foam 
class (at the end of its initialization). The total weight of the event is w{x,i) = Wi{x)/pi. 
Let us call this scenario an "externally organized multibranching" , for short EOM. 

Both methods have certain advantages and disadvantages. In PAID the user does not 
need to organize the optimal/efficient generation of the branching index i. The root cell is 
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divided into N equal-size sub-root cells, which then evolve separately into an independent 
system of cells, adapting individually to the singularities in the i-th component of p. 
Foam adjusts the relative importance of the sub-root cells and their descendants, and 
finds automatically the optimal number of the division cells in the subfoams within the 
requested total memory limit. In the EOM scheme these adjustments for the individual 
branches has to be done by the user. On the other hand, in some cases, the user may 
want to configure the Foam objects for each branch individually. In the EOM scheme, this 
can be done for each Foam object separately. In the PAID scheme this cannot be done, 
because all cells have the same properties, the cell-split algorithm is the same, the cell 
geometry is common, etc. In most cases, the PAID method will be preferred, because it 
is easier to organize. 

In the following we shall concentrate on the PAID scheme. In this case, the normaliza- 
tion integral is provided by the Foam at the end of the exploration phase, and it includes 
the sum over a discrete variable: 

N 

i=i 

We also have the usual relation between the average weight and the integral 

N 

i=l 

The above method extends trivially to the case of several discrete variables. As already 
stressed, the relative probabilities of the discrete components Pi ^ R!^ = J p'^{x)dx^ in the 
MC generation are automatically adjusted by the Foam algorithm, so that the maximum 
weight or the total invariance is minimized. In the EOM scheme, arranging this in the 
user program would require an extra programming effort, while in the PAID method this 
comes for free, by exploiting standard Foam functionality. 



2.11.3 Multilayer method 

There is an alternative PAID-type method of dealing with the problem of the discrete 
variable, which generates points according to p{xi, . . . ,Xn,i), i — 1,2, . . . , N. It will 
produce the same distribution but will differ from PAID in the MC efficiency, in terms of 
the maximum weight or variance. One may simply generate, with the help of Foam, the 
n-dimensional auxiliary distribution 

N 

p{xi,...,Xn) = ^p{xi,...,Xn,j) (9) 

and next, for each generated x, choose randomly a discrete variable i according to the 
probability 

n 

Pi{x) = p{x,i)/^p{x,j). (10) 
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Let us call it PAID*, or the multilayer method. This method is slightly less convenient 
to implement, as is clearly seen for n = 0, where the user effectively has to generate by 
himself the discrete variable i = 1,2, . . . , N according to the above probability, by means of 
creating an inverse cumulative distribution, mapping a random number into i, etc., while 
in the standard PAID scenario all this job is done by the Foam program]^. Furthermore, 
in the PAID method each component distribution p{x,j) may have a "cleaner" structure 
of the singularities than the sum. Consequently, in the PAID method Foam will probably 
find it easier to learn the shape of each component distribution than that of the sum 
in PAID*. These two kinds of equivalent multibranching algorithms, PAID and PAID*, 
are described and analysed in Ref. [^. The PAID* method is used in the KKMC event 
generator of Ref. to generate index i numbering the type (flavour) of the produced 
quark or lepton. 



2.11.4 Multibranching and mapping 



However, the most important reason for setting up Foam according to the PAID scenario, 
with the separate foam build-up for each component distributions p{x,j), is that for each 
component one may apply individually adjusted mapping of variables, which makes ev- 
ery component distribution much less singular. The combination of the mapping and 
multibranching is probably the most powerful known methods of improving MC effi- 
How it can be actually realized with the help of Foam depends on the 



20,21 



ciency 

properties of the distribution p{x) to be generated. In the case where we have an explicit 
sum over many components 



N 



p{Xi, ...,Xn) 



^p(xi. 



each of the components being positive, with a distinctly different and well known structure 
of the singularities, we would recommend the use of Foam in the PAID scheme. Knowing 
the structure of singularities, we may be able to introduce mapping in each component 
separately, which compensates for these singularities with the Jacobian factor. In such a 
case Foam is provided with the following distribution: 



P(l/i, • • • ,1/n, j) = p{x'~(\y), . . .,x\^\y),j) 



dx'^^^ (y) 



dy 



J 



1,2, 



(12) 



understanding that the translation of the discrete index j into a continuous variable yn+i, 
is done in the usual way. The foam of cells is, of course, build-up in the y-variables, 
different for each j-th branch. The user is fully responsible for the proper mapping 
x^^Ky)^ j = 1;2, . . . , A^, and the calculation of the Jacobian factor in every component 
(branch). In the user subprogram providing the p-distribution, the variable yn+i will first 
be translated into index j and then, depending on the value of j, a given type of mapping 
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The mapping Xn+i i is a simple arithmetic operation. 
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will be applied. For the outside part of the code, the index j can be made available, or it 
may be hidden (erased from the record), depending on the needs of a specific application. 

In some cases, however, one does not know any unique way of sphtting the p{x) into 
well defined positive components as Eq. (0), but one may have a rough idea about the 
leading singularities. This means, that one is able to construct the distribution 



N 



Xi, 



(13) 



where p{x,j) have the same type of the leading singularities as p{x), and one controls 
the normalization of singularities in p{x) up to a constant factor; that is for x in the 
neighbourhood of the j-th singular point, a line or a (hyper)plane, only p{x,j) really 
matters, that is p{x) ~ Cjp{x,j), where Cj is not known a priori]^ 

In addition, let us assume that we are able to compensate for the singularities in each 
p{x,j) exactly by dedicated mapping specific to singularities in the j-th branch. In other 
words, the mapping is ideal in each branch: 



dx^^^ (y) 



dy 



P{x,j)' 



(14) 



The above also means, that we know analytically the exact values of the integrals^: 

Rj = f p{x,j)dx"'. 

In such a case we may employ the algorithm of Foam successfully by means of defining 
the "branching ratio" 



= p{yu---,yn,j)/p{yi,---,yn), ^Hy) = 



constructing the distribution to be digested by Foam as 



piyi 



.yn,j) = bj{x)p{x'f\y),...,xli\y)) 



dx^^\y) 



dy 



Y.ip{x{y)^Y 



(15) 



(16) 



and proceeding as in the PAID scheme described previously. 

The role of the function hj{x) is to isolate out from p{x) "a layer" including just one 
known type of singularity. In order to see how this method works, let us consider a 5*^"-' {x — 
a) shape singularity in the j-th component (narrow Gaussian peak etc.) of size e. Then, in 
the vicinity of that singularity we have hj{x) = 1 and p{x) ^ Cjp{x,j), while further away 
we have p{x) ~ (9(e"). The Foam program will, of course, include the Cj factor properly 
in the normalization, and build up the foam of cells everywhere, close to a singularity 
and far away. It will do it really not in the x variables but in the y variables. Now, the 
mapping x ^ y (specific to the j-th branch) will expand the singularity neighbourhood 



^^This definition is not very precise; it roughly means that each component is approximately a product 
of the singular factors and cannot be reduced into the sum of these. 
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In fact, we could normalize p{x,j) to unity, Rj = 1, if we wanted. 
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of 0{e) into a region of size of 0{1), while the y-image of the remaining x-space will be 
of the size 0{e). This can be the source of the following drawback: in the small y-domain 
of 0{e) described above, at the positions of the singularities from the other components 
i 7^ j, p{y,j) may get narrow spikes or dips of height of 0{1), such that their integral 
contribution will be negligible, of 0(e"), and this may not look as a problem. However, 
the Foam algorithm may find it difficult to locate these structures, and this may lead to a 
small but finite bias of the generated distributions and calculated integrals. This should 
be kept in mind, and special tests (MC runs with a maximum number of cells, and high 
MC statistics) should be performed, in order to check that this effect is absent. 

The above method is quite similar to that of Ref. One difference is that in the 
latter there are several iterations with an aim af adjusting the relative normalization of 
the components p{xi, . . . ,Xn,j) to p(x). Our scheme could effectively be regarded as the 
method of Ref. with just one iteration; that is the first step being the foam build-up, 
and the second step (1st iteration) being the MC simulation. One iteration is sufficient in 
the limit of vanishing overlap of the components p{xi, . . . , x„, j) in the entire p{x). While 
in the method of Ref. a better adjustment is provided by the next iterations, in Foam 



the cellular adaptive method provides an extra mileage. One cannot therefore say which 
one is better in general - it depends on the distribution p(x). 

In fact, in the PAID scheme with the mapping, extra iterations are also possible. It 
can be done as follows: (a) read Rj from all "leading cells" after the foam build-up, 
(b) rescale p{yi, . . . ,yn,j) [Rjl Rj)p{.yii ■ ■ ■ yyn,j), and (c) repeat the foam build-up0 
for the new branching ratios bj{x) in Eq. (|1^). The above procedure can be repeated. 
Whether such an iteration is profitable depends on the particular distribution - we expect 
that in most cases it is not necessary, thanks to the adaptive capabilities of Foam. 

Last, but not least, let us also consider the case of a sum of integrals with different 
dimensionality or, in other words, the distribution in which the number rij of the contin- 
uous variables Xi, depends on a certain discrete "master variable" i = I, ...,N (for 
example rii = i): 

N 

R=^ pt{Xi,...,XnJ. (17) 

1=1 

Foam can deal with this case too. The simplest solution is to find the maximum dimen- 
sion njnax and add extra dummy variables on which some of pi does not depend, such 
that formally all sub-distributions have the same dimension nmax- In this way one is back 
to a situation described earlier, and may apply the PAID method, with or without the 
additional mapping. The slight drawback of this solution is that in the present imple- 
mentation of Foam we cannot inhibit the unnecessary cell divisions across the directions 
of the newly introduced dummy variables - simply because they are not the same in all 
branches^. This is the reason why this kind of problem can, in some cases, be dealt with 
more efficiently using the EOM scenario, with a separate Foam object for each branch. 

^^When programming with Foam it is possible to erase all cells from memory and rebuild them. 
^"^A more sophisticated procedure of inhibiting the division could be implemented in Foam, if there is 
a strong demand. 
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For additional practical examples on how to realize multibranching with Foam, see 
Section 



3 Cell split algorithm and geometry 

As already indicated, our algorithm of the cell split covers two strategies: (A) minimization 
of the maximum weight Wmax and (B) minimization of the variance cr, where both Wmax 
and a are calculated in the Monte Carlo generation, using the MC weight w = p/p'. The 
distribution p' is the result of the exploration (it is constant over each cell) and is frozen at 
the end of exploration. The subsequent MC event generation of events is done according 
to p'{x). Its integral R' = J p'{x)d^x has to be known exactly before the start of the 
MC generation. On the other hand, the best estimate of the integral R = J p{x)d"'x is 
obtained, up to a statistical error, at the end of the MC event generation with the usual 
relation: R = R'{w)p'. The average (...)p' is over events generated according to p'. 

There is another important ingredient in the algorithm of the cell split: in addition 
to the auxiliary distributions p'(x) we also define another distribution piossi^) related 
to the integrand p{x). The important role of the distribution piossi^) is to guide the 
build-up of the foam of cells; the function Rioss = J Piossd"x is minimized in the process 
- its value is decreasing step by step, at each cell split. Obviously, both pioss{x) and 
p'{x) are evolving step by step during the foam build-up. Once the division process is 
finished, the distribution pioss{x) is not used anymore; MC events are generated with 
p'{x). Nevertheless, piossix) is strongly related to the properties of the weight distribution 
in the MC generation phase. 

(A) In the case when our ultimate aim is to minimize Wmax, we define 

p'{x) = max p{y), for x G Cellj, 

y^Celli 

f , f 

Rioss = / d"-x [p (x) - p{x)] = / d'^x pioss{x). 

The distribution pioss is the difference between the "ceiling distribution" p' and the actual 
distribution p from which it is derived. The rejection rate in the final MC run is just 
proportional to the integral over the loss distribution piossi^) by construction, i.e. the 
rejection rate = Rioss/ R- (This justifies the name "loss".) The distribution pioss{x) also 
has a clear geometrical meaning, see below. 

(B) In the case when we do not care so much about the maximum weight and the 
rejection rate, but we want rather to minimize the ratio of the variance to the average of 
the weight, cr/ (w), in the final MC generation, we are led to the following definition: 

p'i^) = \fW)i^ for X e Celli, 

I (l^J 

Pioss{x) = V (p2)^ - (p)^, for s G Celli. 

The average (...)/ is over the J-th cell; see Appendix A for a detailed derivation of the 
above prescription. The ratio a / (w) in the final MC generation is a monotonous ascending 
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function of Rioss = J Pioss{x)dx"', see Appendix A. Consequently, minimization of Rioss is 
equivalent to minimization of a/ (w). 

3.1 Rules governing the binary split of a cell 

The basic rules governing the development of the foam of cells are the following: 

(a) For the next cell to be split we chose a cell with the largest^ Rioss- 

(b) The position/direction of a plane dividing a parent cell into two daughter cells u — *• 
u' + oj" is chosen to get the largest possible decrease ARioss = Rioss ~ ^toss ~ ^toss- 
How is the split of a given cell into two daughter cells in step (b) done in practice? The 
method relies upon a small MC exercise within a cell, in which a few hundreds of events 
are generated with a flat distribution. They are weighted with p{x) and projected onto 
k (hyperrectangular case) or n{n + l)/2 (simplical case) edges of the cells. In the mixed 
case of a cell being the Cartesian product of a fc-dimensional hyperrectangle and an n- 
dimensional simplex, there are k + n{n + l)/2 projections/edges. Resulting histograms 
are analysed and the best "division edge" and "division hyperplane position" are found - 
the one for which the estimate (forecast) of the ARioss is the largest. In the actual Foam 
algorithm, each new-born cell is immediately explored, its Rioss, R and R' are calculated, 
and the best candidate as to the direction and position of the dividing plane are established 
and memorized, as the attributes of the cell; see below for details. In this way, every newly 
created cell is ready for an immediate binary division. 

3.2 Geometry of the binary split of a cell 




Figure 5: Geometry of the split of a S-dimensional cell, simplex or hyperrectangle. 

^^In the Foam code there is also an option of choosing the next ccU to be spht randomly, according to 
probability proportional to Rioss, instead of a cell with the largest Rioss- 
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In Fig. 1^ we show a 3- dimensional cell being a simplex or a hyperrectangle and we 
visualize the geometry of their split. 

Let us first describe the split of the n-dimensional simplical parent cell into two daugh- 
ter cells. In this case [0 we need to know which of n{n + l)/2 edges defined by any pair 
of vertices of a given simplex is the best for the split. Suppose that it is an edge defined 
by a pair of indices i 7^ j, where i,j = 1,2, Ny, of the vertices (yKi,VKj), see 



Sect. ^]4|for the method of numbering the vertices. A new vertex Vny+i is put somewhere 



on the line (edge) in between the two vertices: 

V^y+i = XVk, + {1 - \div)VK„ 0<Xdiv<l, (20) 

where the division parameter X^iv is determined by using an elaborate procedure described 
later in this section, and the number of vertices is updated Ny —>■ Ny + 1. With the new 
vertex two daughter simplices are formed with the following two lists of vertices (their 
pointers): 

(iTi, K2, K,_^, {Nv + 1), Ki+^, K,^i, Kj, Kj+^, Kn, K^+i), 
(ifi, K2, Ki.i, K„ Ki+i, (Nv + 1), Kj+i, Kn, Kn+i). ^ ' 

For the fc-dimensional hyperrectangular cell defined with a pair of vectors {q, h) we 
first decide on the direction of the division hyperplane. Assuming that this hyperplane 
is perpendicular to the i-th direction, the two daughter cells (a) and (b) are defined with 
the two pairs of new vectors as follows: 

= (91, 92, • • • , qk), h!-"^ = {hi, h2,..., hi-i, hiXdiv, hi+i, ...,hk), 

<t~^^ = {Qij • • • 5 Qi-ij Qi + hiXdiv, hi+i, . . . ,qk), hS^^ = {hi, . . . , /ij-i, hi{l — Xdiv), /^j+i, • • • , h^) 

(22) 

The 3-dimensional case of the simplical and hyperrectangular cell split made in this 
way is illustrated in Fig. |. 

3.3 Projecting points into an edge 

Before we describe the determination of the division edge and of the division parameter 
Xdiv, let us still discuss certain geometric aspects of the Foam algorithm - that is how we 
project a point x inside a cell onto one of the edges of the cell. In the case of a simplex 
the edges are numbered by the pair of indices {i,j), i > j, which number edges spanned 
by a pair of vertices^ {Vk^^Vkj), while in the case of the hyperrectangle the i-th. edge 
is spanned by the pair of vectors q and g+ = {qi, . . . , qi-i, qi + hi, . . . , qn). The point x 
inside a cell is projected onto the edge and parametrized using the parameter A G (0, 1). 
The parameter A will be used to define an auxiliary projection dp/dX for each edge in the 



^^As explained in Sect. 2.4, the numbering of vertices is done using pointers Ki to the elements of the 
array of vertices. 
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(1,2) 



Figure 6: Geometry of the split of the 3-dimensional simplex eel 





Figure 7: Geometry of the split of the 3-dimensional hyperrectangular eel 



following subsection. In particular we have to know how to evaluate A in an efficient way. 
In the case of a simplex cell we have: 



pro] 



XVK^ + {l-\j)VK,, 0< Ai,- < 1, 



(23) 



where 



Xij{x) 





Det, 




|Deti 


1 + 1 


Detjl 



(24) 



Deti = Det(fi, . . . , r;_i, fj+i, . . . , r;, r^+i), 

Detj = Det(fi, . . . , rj_i, fj+i, . . . , r„, f„+i), = Vk^ - x, 

and Det(xi, X2, ■■■,Xn) is determinant. The case of a hyperrectangular cell is much simpler: 

Xi = {xi - qi)/hi. (25) 

Obviously, owing to the time-consuming evaluation of the determinants at higher dimen- 
sions, the above projection procedure will be much slower for simplices than for hyper- 
rectangles. 
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In Fig. 1^ we illustrate a projection procedure into six edges for 3-dimensional simplex 
and in Fig. ^ the case of the three edges of the 3-dimensional hyperrectangular cell. 
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Figure 8: Projection histogram. The case of optimizing the maximum weight. 



3.4 Determination of an optimal division edge and of X^iv 

Our aim is to find out which division hyperplane, cutting through which edge, provides 
the best gain of the total integral Rioss, summed over two daughters, as compared to the 
parent cell. In order to do that, we first analyse all possible division planes, for all edges, 
and find out the best one, in terms of the gain in Rioss- In other words, we go through 
all edges {k edges for a hyperrectangle and/or n{n + l)/2 for a simplex) and find, for 
each edge, the optimal parameter X^iv and the corresponding best gain in Rioss- We then 
compare between the gains in Rioss for all edges and define the optimal edge as the one 
with the best gain in Rioss- The procedure of finding the best Xdiv is essentially the same 
for simplical and hyperrectangular cells; on the other hand, in the algorithm for finding 
the best X^iv there is a difference between the cases of optimization of the maximum 
weight and of the variance; see the following discussion. 



3.4.1 Optimization of the maximum weight — choosing A 



div 



Let us consider first the case of finding the best Xdiv for Rioss corresponding to the opti- 
mization of the maximum weight. The 1-dimensional case is a good starting point. The 
cell in this case is just an interval {q,q + h) and X = {x — q) /h. In the left-hand part 
of Fig. ^, we see a histogram with Ni, bins, made of 1000 events generated inside a cell 
(interval) using the weight w = p{x), so that the histogram represents approximately the 
distribution dp/dX, X G (0,1). This distribution (histogram) peaks close to the lower 
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edge. The function p\x) = max^^gcez/ p{x) is constant over the cell and is depicted as an 
upper horizontal line marked "Parent p'" . The contribution of this particular (parent) 
cell to Rioss = J^^iiPiossdx = J^^iiip'ix) — p{x))dx is easily recognized as an area between 
the line marked "Parent p'" and the histogram line. If we have stopped the exploration 
at this stage, with this parent cell, then in the MC run points would be generated with 
the flat "Parent p'" and the weight would he w = p{x) / p. Turning weighted events into 
unweighted ones by accepting r < w events and rejecting r > w, where < r < 1 is a 
uniform random number, would correspond to generating points (A,r) within the rectan- 
gle below the "Parent p'" line, accepting all points below the histogram line and rejecting 
all those above, in the area marked "Parent Rioss" ■ This justifies the subscript "loss". 

The best cell division is found by examining all A*";, — 1 end-points A = q + ih/Nf,, 
i = 1,2, Nb — 1 of the bins in the histogram, as a possible candidate for the division 
point (plane in two and more dimensions) between the two daughter cells, and choosing 
the best one. In the right-hand part of Fig. |^ we have marked such a candidate division 
point with a star. For a given division point, we determine for two daughter cells the new 
"ceiling function" p'; in Fig. ^ it is the line marked "New p'" . For each daughter cell we 
evaluate Rioss- The summary Rioss for both daughter cells is easily recognized as an area 
between the line marked "New p'" and the histogram. Of course, we automatically get 
the new total Rioss smaller than for the original parent cell! This procedure is repeated 
for all possible j = 1,2, Nf, — 1 division points and each time we record the net gain in 

^jRloss = Rloss,parent — Rioss, daughterl — Rioss, daughter2- For the actual best divisioU poiut We 

chose the division point with the largest gain A j Rioss- In Fig. || the star marks the best 
division point. 




Figure 9: Two-dimensional p{x) (left) and the geometry of the first three simplical cells (right). 
Inside the area marked by the inset rectangle p{x) = 0. 
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Figure 10: Distributions used in the construction of Xdiv in the case of optimizing the maximum 
weight. 



Let us now consider the 2- dimensional distribution p{x) depicted in the left part of 
Fig. ^, which is non-zero within the narrow strip along the four edges of the rectangle. In 
the simplical mode Foam divides starting ra-dimensional hyperrectangle into n\ simplices 

- in this case into two triangles. We concentrate on the division procedure of the lower 
triangle; see the right part of Fig. |^. In the exploration of this triangular cell we use 1000 
MC points and project them onto three edges. The corresponding three histograms are 
shown in Fig. |10|, where the middle histogram represents the projection onto the diagonal 

- this is why it features two distinct peaks at the ends. In all three plots we have also 
drawn the curve for p'{x) for the best hypothetical split. The most promising split in 
terms of the gain in Rioss turn outs to be related to the middle plot and is marked in the 
right part of Fig. |^. 

The reader may notice that the p'{x) in the middle plot of Fig. ^ is not of the type 
discussed above, because it has two discontinuities instead of one. This is because in Foam 
we have introduced an important refinement of the algorithm to find an optimal A^^. 
One may easily notice that the algorithm described above could not correctly locate the 
drop in the distribution dp/dX of the middle plot, because there are two equally strong 
peaks at the ends of this distribution^. In the improved version of the algorithm the 
search of the optimal Xdiv uses all pairs of the bin edges (Aj, Xj) = {q + ih/Nb, q + jh/Nb), 
< < J < Nfy. For every pair a new ceiling function p'{x) is determined, such 
that it is unchanged outside the subinterval (Aj, A^) and is "majorizing" the histogram 
bins inside this subinterval. Once we find out the best pair in terms of Rioss, then 
we take either Aj or Aj as a division point X^iv (at least one of them is not equal to or 
1). In the case of two or more peaks in dp/dX the resulting division point A^j^ happens to 
be close to one of the edges of the gap between the two peaks. This feature prevents the 
Foam algorithm (at least partly) from placing a new division plane across a void in the 
multidimensional distribution p{x). In other words such a void will "repel" the division 
hyperplanes. In the case of the double peak structure in the middle plot of Fig. the 
improved algorithm will of course allocate the big value of Rioss to a new cell (interval). 
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The algorithm would pick up Xdiv at a random point between the two peaks. 
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which includes the gap. In the next step of the foam build-up this cell (interval) will have 
a big chance to be split, and for this split the position of the split point will be located at 
the second edge of the gap. This is exactly what we need for an efficient foam evolution. 



of 
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Figure 11: Distributions used in the construction of A^j^ in the case of optimizing the variance. 



3.4.2 Optimization of the variance — choosing A 



div 



Let us consider now the case of finding out the best X^iv for Rioss corresponding to op- 
timization of the variance. The strategy is again to choose Xdiv by minimizing Rioss- In 



Fig. |IT] we illustrate our algorithm on the example of the three projections of the trian- 
gular cell. The three projections correspond to a triangular cell in two dimensions. (We 
do not specify p{x), as it is irrelevant for the purpose of our explanation.) In the upper 
row of the three plots in Fig. ^ we show as a horizontal line the value of the distribution 
Pioss = a/ (p^) (it is the same for all three projections). The solid histogram is the distri- 
bution dp/d\ and the dashed histogram is the distribution dpioss/dX calculated bin by bin 
using treating every bin as a separate cell. The properly normalized difference of 

the above two distributions — (p) is plotted separately as the dashed histograms in 

the lower row of the three plots in Fig. |ll|. The horizontal line for the total pioss is shown 
once again there. The histogram of dpioss/dX gives us an idea, where the biggest source 
of the variance is located and our aim is to "trap" it properly with an intelligent choice of 
Xdiv We follow an algorithm similar to that of the maximum weight minimization, namely 
we loop over pairs of the bin edges (Aj, Xj) = {q + ih/Ni,, q + jh/Nf,), < i < j < N^, 
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and for every pair we calculate Rioss inside the interval (Aj,Aj) and outside it. We find 

out which provides the biggest gain AijRioss = Rloss,parent- Rloss,Inside- Rloss,Outside- 

In the lower row of the plots in Fig. 11 we show as a solid line the distribution of Rioss 
for the best pair Depending on the peak structure, at least one of the division 

points of the optimal pair (Aj, Xj) is different from or 1, and we take this one as Xdiv 
In Fig. |Tl] the chosen X^m are marked with black triangles. The above procedure is done 
for each edge, and the best AijRioss is used as a guide to define an edge for which the 
next cell division will be executed. The information about the best edge and the best 



division point Xdiv is recorded in the cell object. As seen in Fig. 11, Xdiv tends to fall at 



the position where dpioss/dX drops or increases sharply. Note that since the division point 
is always at the edge of the bin, it is therefore a rational number, Xdiv = j/Nt. This has 
interesting consequences, since the number of the bins iVf, is fixed, it is therefore sufficient 
to memorize this integer index j (2 Bytes) together with the integer index of the division 
cell edge (also 2 Bytes) as an attribute of the cell, in order to define fully and uniquely 
the geometry of the division of the cell! See Section for more details on how this is 
exploited to save the computer memory needed to encode the entire foam of cells. 



3.4.3 Concluding remarks on the cell division algorithm 

The algorithm of the split of the cell is the important and most sophisticated part of the 
new Foam. Let us therefore add a couple of final remarks: 

• The new, much improved procedure of the choice of the division plane is the most 
significant difference[^ with respect to Foam 1.x of Ref. [|1|. 

• The choice of the edge based on the histograms for each edge makes sense if we use 
histograms with at least 2-4 bins and at least 100 MC events per cell. This might 
be a serious limitation for these p{x) , which require a lot of CPU time per function 
call. 

Finally, in Fig. |1^ we show examples of the evolution of the foam of cells as the number 
of the cells grows gradually. The 2-dimensional case is easily visualized and we show it in 
Fig. for triangular and rectangular cells. In the upper six plots, p(x) features a circular 
ridge, in the bottom two plots it is concentrated along the antidiagonal Xi + X2 = 1, the 
last one corresponding to the p{x) of Fig ||. 



3.5 Limitations 

We are aware that the present procedure of selecting the next cell for the split and the 
procedure of defining the division plane, although quite sophisticated, is not a perfect one 
and has certain shortcomings. Some of them can probably be corrected for, but some of 
them are inherent to the method. In Fig. |13] we show a surprisingly simple example of a 
function for which our method of finding a good division point Xdiv fails. It fails simply 



I would like to thank A. Para for a discussion which ignited this new development. 
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Figure 12: Examples of the 2-dimensional foam. Number of cells from 10 to 2500. 



because both projections of p{x) onto two edges of the rectangle are just flat and our 
procedure will pick up some A^j^, randomly within (0, 1), while the most economic division 
point is in the middle X^iv = 1/2- On the other hand, although the Foam algorithm gets 
disoriented for the first division, it will recover and correct for the false start in the next 
divisions rather quickly. It will eliminate the two voids from its area of interest. 

Let us notice that the distribution of Fig. |TB| violates maximally the "principle of 
factorizability" p[x) = YYi Pi{^i)i the principle on which the VEGAS family of the programs 
is built Contrary to VEGAS the problem with factorizability in Foam is not a general 

one, but is limited to a single cell and usually goes away after the cell split. Nevertheless, 
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Figure 13: Example of p{x) for which the Foam algorithm of cell division fails. 



the algorithm of Foam analysing projections on all edges in a single cell is relying on the 
"principle of factorizability" . 



Class 


Short description 


TFDAM_INTEGRAND 
TFVECT 

TFMATRIX 

TFPARTITION 

TFCELL 

TFOAM 


Abstract class (interface) for the Foam integrand distribution p(x) 
Class of vectors with dynamic allocation of the components. Used in 
TFOAM and TFCELL 

Square matrices, used for simplical geometry in the foam 

Auxiliary small class for looping over partitions and permutations 

Class of objects presenting single cell used in TFOAM (Cartesian product 

of the simplex and hyperrectangle) 

Main class of Foam. The entire MC simulator 


TPSEMAR 
TFHST 

TFMAXWT 
TFDISTR 


Marsaglia etal. random number generator 23]. 

Simple class of one-dimensional histograms. Used only in the Foam ver- 
sion without ROOT 

Monitors MC weight, measures performance of the MC run 
Collections of distributions p{x) for testing Foam 



Table 1: Description of C++ classes of Foam. 



4 The Foam code 

At present, the C++ version of the Foam code is more advanced than the Fortran?? 
version. (We do not plan to develop F77 code any further.) In this section we shall 
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TFOAM member 


Short description 


float m_Version^ 
char m_Date[40] 
char m_Name[128] 
int m_nDim** 
int m_kDim* 
int m_TotDim^ 
int m_nCells'* 
int m_vMax 
int m_LastVe 
int m_RNmax 


Actual version of the Foam (like 2.34) 

Release date of the Foam 

Name of a given instance of the TFOAM class 

Dimension of the simplical subspace 

Dimension of the hyperrectangular subspace 

Total dimension = m_nDim+m_kDim 

Maximum number of cells 

Maximum number of vertices (calculated) 

Actual index of the last vertex 

Maximum number of random numbers generated at once 


int m_OptDrive* 
int m_OptEdge* 
int m_OptPeek* 
int m_OptOrd'* 
int m_OptMCelP 
int m.Chat** 
int m_OptDebug 
int m_OptCulst 
int m_OptRej 


Type of optimization =1,2 for variance or maximum weight reduction 

Decides whether vertices are included in the cell MC exploration 

Type of cell peek =0,1,2 for maximum, random, random2 

Root cell is simplex for OptOrd=l, hyperrectangle for OptOrd=2 

=1 economic memory for hyper rectangles is on; =0 off 

=0,1,2 chat level in output; =1 for normal output 

=1, additional histogram (dip-switch) 

=1, numbering starts with hyperrectangle (dip-switch) 

=0 for weighted events; =1 for unweighted events in MC generation 


int m_nBin'* 
int m_nSampF 
int m_EvPerBin'^ 
int m_nProj 


No. of bins in edge histogram for cell MC exploration 
No. of MC events, when dividing (exploring) cell 
Maximum number of effective {w = 1) events per bin 
Number of projection edges (calculated) 



Table 2: Data members of the TFOAM class. Associated setters and getters marked as super- 
scripts s and g. 



therefore describe mainly the C++ code. 

The code of the Foam version 1.x was originally written in Fortran?? with popular 
language extensions, such as long variable names etc. It was already written in an object- 
oriented style, as much as it had been possible. In particular the main classes TFOAM 
and TFCELL of the present C++ version were already present under a certain form. The 
important shortcoming of the F77 version is the lack of dynamic allocation of the memory. 
Otherwise, it has most of the functionality of the C++ version; see later in this section 
for a list of limitations. 



4.1 Description of C-\ — |- classes 

In Table |1] we list all classes of the Foam package. The main class is TFOAM, which is the 
MC simulator itself. It is served by the class TFCELL of the cell objects, and three auxiliary 
classes TFVECT, TFMATRIX and TFPARTITION. The other classes are not related directly to 
the Foam algorithm - they are utilities used by Foam: random-number generator class 
TPSEMAR [^] and the histogramming class TFHIST. The class TFDIST provides a menu of 



33 



TFOAM member 


Short description 


Provision for llic; mulliV^raiicliiiig 


int *m_MaskDiv 
int *m_InhiDiv 
int m.OptPRD 
TFVEC^r **m_XdivPRD 


![m_nProj] Dynamic mask for cell division 
![ni_kDim] Flags inhibiting cell division, h-rcctang. subspace 
Option switch for predefined division, for quick check 
ILisIs of division values cncodcxl in one vector p(;r direction 


Geometry of cells 


int m_NoAct 
int m_LastCe 
TFCELL **m_Cells 
TFVECT **m_VerX 


Number of active cells 

Index of the last cell 

[m_nCells] Array of ALL cells 

[m_vMax] Array of pointers to vertex vectors 


Monte Carlo generation 


double m_MaxWtRej; 
TFMAXWT *m_MCMonit; 
TFCELL **m_CellsAct 

double *m_PrimAcu 
TObjArray *m_HistEdg 
TObjArray *m_HistDbg 
THID *m_HistWt; 
TFHST **m_HistEdg 
TFHST *m_HistWt; 


Maximum weight in rejection for getting w = 1 events 
Monitor of the MC weight for measuring MC efficiency 
!Array of pointers to active cells, constructed at the end of 
foam build-up 

! Array of cumulative Yli=i cell index generation 
Histograms of w, one for each edge, with ROOT 
Histograms for debug (m_OptDebug=l), with ROOT 
Histograms of MC weight, with ROOT 
Array of pointers to histograms, without ROOT 
Histograms of MC weight, without ROOT 


double *m_MCvects 

double m_MCwt^ 
double *m_Rvec 


[m_TotDim] Generated MC vector for the outside user 

MC weight 

[m_RNmax] Random number vector from r.n. generator, up 
to m.TotDim+l maximum elements 


Exloruals 


TFOAMJNTEGRAND *m_Rho^s 
TPSEMAR *m_PseRan5s 


The distribution p to be generated/integrated 
Generator of the uniform pseudo-random numbers 


Statistics and MC results 


long m_nCalls^ 

long m_nEfFev 

double m_SumWt, m_SumWt2 

double m_NevGen 

double m.WtMax, m.WtMin 

double m_Prinic^ 

double m_MCresult 

double m_MCerror 


Number of function calls 

Total no. of effective w = 1 events in build-up 
Sum of weight w and squares 
No. of MC events 

Maximum/Minimum weight (absolute) 

Primary integral R', (i? = R{w)) 

True integral R from the cell exploration MC 

and its error 


Worl 


dng space for cell exploration 


double *m_Lambda 
double *m_Alpha 


[m_nDim] Internal parameters of the simplex: ^ < 1 
[m_kDim] Internal parameters of the h-rectang.: < CKj < 1 


Table 3: Data members of the class TFOAM. Cont. 
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TFOAM method 


Short description 


Constructors and destructors 


TFOAM 

TFOAM(const char*) 
TFOAMQ 

TFOAM(const TFOAM&) 

TFOAM& operator = (const TFOAM& ) 


Default constructor (for ROOT streamer) 

User constructor 

Explicit destructor 

Copy Constructor NOT USED 

Substitution NOT USED 


Initiahzation, foam build-up 


void Initiahze(TPSEMAR*, 

TFOAMJNTEGRAND*) 

void Init Vertices (void) 

void InitCclls(void) 

void Grow(void) 

int Divide(TFCELL *) 

void Explore(TFCELL *Cell) 

void Carver (int&,double&,double&:) 

void Varedu(double[ ],int&;,double&;, doubled) 

long PeekMax(void) 

TFCELL* PeekRan(void) 

void MakcLambda(void) 

void MakeAlpha(void) 

int CeUFill(int, TFCELL*, 

int*,TFVECT*,TFVECT*) 

void MakeActiveList(void) 


Initialization, allocation of memory 
and the foam build up 
Initializes first vertices of the root cell 
Initializes first n! cells in h-rect. root cell 
Adds new cells to foam, until buffer is full 
Divides cell into two daughters 
MC exploration of cell main subprogram 
Determines the best edge, liJmax-reduction 
Determines the best edge, cr-reduction 
Chooses one active cell, used in Grow 
Chooses randomly one active cell, in Grow 
Generates random point inside simplex 
Generates rand, point inside h-rectangle 

Fills next cell and return its index 
Creates table of all active cells 


Generation 


void MakeEvent(void) 
void GetMCvect (double *) 

void GetMCwt (double &) 
double MCgenerate(doublc *MCvect) 
void GenerCell(TFCELL *&) 
void GenerCel2 (TFCELL *&) 


Makes (generates) single MC event 
Provides generated random MC vector 
Provides MC weight 
All the above in single method 
Chooses one cell with probability ~ R'^ 
Chooses one cell with probability ~ R'^ 


Finalization, reinitialization 


void Finalize(doublc&:, doublcfe) 
void GetlntegMC (doubled, doubled) 
void GetIntNorm(double&;, doubled) 
void Get WtParams (const double, 

doublccfe, doublefc, doubled) 
void LinkCells(void) 


Prints summary of MC integration 
Provides MC integral 
Provides normalization 

Provides MC weight parameters 
Restores pointers after restoring from disk 


Debug 


void Check All (const int) 
void PrintCells(void) 
void Print Vertices (void) 
void LaTexPlot2dim(char*) 
void RootPlot2dim(char*) 


Checks correctness of the data structure 
Prints all cells 
Prints all vertices 

Makes LaTeX file for drawing 2-dim. foam 
Makes C-|— I- code for drawing 2-dim. foam 



Table 4: Methods of TFOAM class. 
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the distributions for testing Foam. In the following we shall describe in more detail the 
key classes TFOAM and TFCELL. 

4.2 TFOAM class 

TFOAM is the main class. Every new instance of this class (properly initialized) is another 
independent Foam event generator. In Tables ^ and ^ we provide a full list of the data 
members of the class TFOAM and their short description. As seen in these tables, we have 
added the prefix "m_" to all names of the data members, so that they are visually different, 
in the code, from the other variables (which is recommended practice in C++ coding). 

Generally, one may notice that many data members could be declared (and allocated) 
as the local variables in the procedures, instead of being data members. For example, 
vector m_MCvect, which transpors random numbers out from the random number generator 
m_PseMar could be declared locally at every place where m_PseMar is called. We opted for 
a more "static" structure of the data, with more than necessary of the data members in 
the class, at the expense of the human readability of the code, in order to: (a) facilitate 
the implementation persistency with ROOT, (b) gain in execution speed and (c) facilitate 
the translation to other languages. 

Most of the methods (procedures) of the class TFOAM are listed in Table ^ We omitted 
in this table "setters" and "getters", which provide access to some data members, and 
simple inline functions, such as sqr for squaring a double variable. Data members which 
are served by the setters and getters are marked in Tables and | by the superscripts 
"s" or/and "c/". 

Let us now characterize briefly the role of most important methods of the class TFOAM 
in the Foam algorithm. 

4.2.1 Procedures for Foam initialization and foam build-up 

The constructor TFOAM (const char*) is for creating an object of the class TFOAM. Its 
parameter is the name given by the user to an object. The principal role of this constructor 
is to initialize data members to their default values - no memory allocation is done 
at this stage. After resetting all kinds of steering parameters of the Foam to preferred 
values (using setters) the user calls the Initialize method, which builds up the foam 
of cells. The two methods InitVertices and InitCells allocate arrays of vertices and 
cells (pointers) with empty cells. The empty cells are allocated/filled using CellFill. 
Next comes the procedure Grow which loops over cells, picking up the most promising 
cell for the split, either randomly using PeekRand or deterministically using Peekmax. 
The chosen cell is split using Divide. It is, however, the procedure Explore called by 
Divide (and by InitCells for the root cell) which does the most important job in the 
Foam build-up: it performs a small MC run for each newly allocated daughter cell. It 
calculates how profitable the future split of the cell will be and defines the optimal cell 
division geometry with the help of Carver or Varedu procedures, for maximum weight 
or variance optimization respectively. All essential results of the exploration are written 
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InitVertices 




InitCells 




Grow 



PeekMax 




PeekRand 




Divide 



CellFill 



I Explore | 



Figure 14: Calling sequence of the Foam procedures during the foam build-up (initialization). 

into the explored cell object. At the very end of the foam build-up, MakeActiveList is 
invoked to create a list of pointers to all active cells, for the purpose of the quick access 
during the MC generation. The procedure Explore uses two procedures MakeLambda and 
MakeAlpha, which generate randomly (uniformly) coordinates of the MC points inside a 



given cell. The above sequence of the procedure calls is depicted in Fig. [14. 
4.2.2 Procedures for MC generation 

The MC generation of a single MC event is done by invoking MakeEvent, which chooses 
randomly a cell with the help of procedure]^ GenerCell2 and, next, the internal coor- 
dinates of the point within the cell using MakeLambda and/or MakeAlpha. The absolute 
coordinates of the MC event are calculated and stored in the data member double-precision 
vector m_MCvect. The MC weight is calculated using an external procedure that provides 
the density distribution p(x), which is represented by the pointer m_Rho. A class to which 
the object m_Rho belongs must inherit from the abstract class TF0AM_1NTEGRAND. The MC 
event (double-precision vector) and its weight are available through getters GetMCvect 
and GetMCwt. Note that the variables of the hyperrectangular subspace come first in the 
m_MCvect, before variables of the simplical subspace. 

The user may alternatively call MCgenerate, which invokes MakeEvent and provides 
a MC event and its weight, all at the same time. 



■^''^Method GenerCelll exists, but it is not used. 
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4.2.3 Procedures for finalization and debugging 

The use of the method Finalize is not mandatory. It prints statistics and calculates 
the estimate of the integral using the average weight from the MC run. The amount of 
printed information depends on the values of m_chat. For the normalization of the plots 
and integrals, the user needs to know the exact value of R' = J p'{x)dx, which is provided 
by the method GetlntNorm or Finalize. The actual value of the integrand from the MC 
series is provided by GetlntegMC. Note that for the convenience of the user GetlntNorm 
provides R' or MC estimate oi R = j p{x)dx, depending on whether the MC run was 
with weighted or w = 1 events. 

Another useful finalization frocedure, 

GetWtParams (const double eps, double &AveWt, double feWtMax, double feSigma) 

provides three parameters that characterize the MC weight distribution: the average 
weight AveWt, the "intelligent" maximum weight WtMax = w^^^ for a given value of 
eps = e (see Sect. |] for its definition) and the variance sigma = a. In particular, in the 
case of unweighted events, if^ax '^^^ ^e used as an input for the next MC run. 

The Foam program includes procedure CheckAll, which checks the correctness of the 
pointers in the doubly linked tree of cells (this can take time for large Nc). It can 
sometimes be useful for debugging purposes. Another two methods PrintVertices and 
PrintCells can be used at any stage of the calculation in order to print the list of all 
cells and vertices. In the case of the two-dimensions there is a possibility to view the 
geometry of the cells with a 2-dimensional plot, which is either a LaTeX file produced by 
LaTexPlot2dim, or a ROOT file produced by RootPlot2dim. 



4.3 TFCELL class 

TFCELL is an important class of objects representing a single cell of the foam. Data 
members of the class are listed in Table ^. 

Most of the methods of the TFCELL class are setters and getters. The non-trivial 
methods are GetHcub and GetHSize, which calculate the absolute position and size of 
hyperrectangles in the algorithm of Section p.6| , and MakeVolume, which calculates the 
Cartesian volume of the cell. In the simplical subspace, the volume is a determinant of a 
square matrix of the class TFMATRIX. 



4.4 Persistency with the help of ROOT 

The C++ language does not provide any built-in mechanism of the persistency of the 



classes. For this purpose we use the ROOT package |T6[, with the help of its "automatic 



streamers". ROOT is a useful C++ library for histogramming, organizing large databases 
of identical objects of the type used in high energy physics experiments. It also provides 
an efficient input /output, with compressing capabilities. 

Providing full persistency for C++ classes preserving all the structures of the pointers 
is not a trivial task in general. ROOT can do it, even for pointers. For this, our code had 
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TCELL member 


Short description 


"Static" members, the same for all cells! 


short int m_kDim 


Dimension of hyperrectangular subspace 


short int m_nDim 


Dimension of simplical subspace 


short int m_OptMCell 


Option of economic memory for usage (hyperrectangular sub- 
space) 


short int m_OptCulst 


=1, Numbering of dims starts with hyperrectangles; =0 simplices 


int m_nVert 


No. of vertices in the simplex = m_nDim-|-l 


TFCELL **m_CellO 


! Pointer of the root cell 


TFVECT **m_VertO 


! Pointer of the vertex list 


Linked tree organization 


int m_Serial 


Serial number (index in m_CellO) 


int m_Status 


Status (active or inactive) 


int m_Parent 


Pointer to parent cell 


int m_DaughtO 


Pointer to daughter 1 


int m_D aught 1 


Pointer to daughter 2 


The best split geometry from the MC exploration 


double m_Xdiv 


Factor A of the cell split 


int m_Best 


The best edge candidate for the cell split 


Integrals of all kinds 


double m_Volume 


Cartesian volume of this cell 


double mJntegral 


Integral over cell (estimate from exploration) 


double m_Drive 


Driver integral Rioss for cell build-up 


double m_Primary 


Primary integral R' for MC generation 


Geometry of the cell 


int *m_ Verts 


[m_nVert] Pointer to array of vertices in simplical subspace 


TFVECT *m_Posi 


Pointer to position vector, hyperrectangular subspace 


TFVECT *m_Size 


Pointer to size vector, hyperrectangular subspace 



Table 5: Data members of the class TFCELL. 



to be reorganized slightly: puting additional directives for ROOT in the comment lines, 
remooving static variable, explicit integer indices instead of pointers in some places. As 
a whole, this solution is relatively simple, and works correctly. In Tables ^ and ^ we have 
at the beginning of the description certain characteristic marks which are directives for 
the persistency mechanism of ROOT, see the ROOT manual for more details. 



One has to remember, when reading a TFOAM class object from the disk, that the 
method LinkCellsO has to be invoked in order to fully reconstruct all pointers in the 
doubly linked tree of cells. Moreover, any object of the class TFOAM restored from the 
disk file will have its internal object for the random number generator and distribution 
function. There is a method that provides access (pointer) to these objects, if necessary. 
The relevant fragment of the code may look as follows: 

TPSEMAR *RNGen= FoainX->GetPseRan() ; //get pointer of RN generator 

TFDISTR *RHO = (TFDISTR*)FoainX->GetRho() ; //get pointer of distribution 
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It might be useful if, for instance, one wants to reinitialize the random number generator 
used by the TFOAM class object, which has been read from the disk file, with a new "random 
seed". 

Foam can be used with or without ROOT. In the code all parts of the code depen- 
dent on ROOT enclosed in the pair of preprocessor commands #if def ROOT_DEF . . . 
#endif , where the ROOT_DEF environmental variable is defined centrally in the header file 
ROOT_DEF.h. Eliminating ROOT requires removing this variable and modifying makefile 
accordingly (the TFHST class has to be linked). The version without ROOT does not fea- 
ture persistency, and is employing its own simple histogramming class TFHST instead of 
the ROOT class THID. ROOT also helps to create documentation of the Foam in the html 
format. We recommend a version tied up with ROOT. 



4.5 Fortran77 version and its limitations 

We also provide users with the Fortran?? versions of Foam. There are two of them at 
the development level 2.02 (May 2001) of the algorithm. In the first one cells can be 
simplical, hyperrectangular and the Cartesian product of the two. This version is limited 
to dimension five for the simplical subspace and not very useful for large dimensions 
(n > 5) in the hyperrectangular space, because it does not feature the memory-saving 
algorithm of Section [2^ . Another version, named MCell (standing for Mega-Cell), features 
only hyperrectangular cells; on the other hand, it includes the memory-saving algorithm 



described in Section p.6| . We recommend the reader to use the version MCell. 

Neither of these versions can have dynamic memory allocation; they have a maxi- 
mum dimensions of the integration/simulation subspaces (simplical and hyperrectangu- 
lar) hard-coded in the source code. Any change of these maximum dimensions requires 
recompilation of the code. 

Present versions in Fortran?? are substantially improved with respect to the original 
version of Ref. For the option of minimizing the maximum weight, they have exactly 
the same algorithm (of the cell split) as the C++ version. They feature, however, an 
older, more primitive version of the algorithm of finding the best cell division for the 
variance reduction. 

The structure of the programs, naming of procedures and variables, configuration 
parameters and their meaning are very similar in the F?? and C++ versions. Some 
differences in the usage will be indicated in the next subsections. 



4.6 Future development 

In the following we indicate some of the possible future developments of the Foam package. 
As already indicated we do not plan to develop the Fortran?? version any further. On 
the other hand, it would be interesting to upgrade the existing Foam version 1.x in the 
JAVA language to the level of the present version 2.x. 

As for the C++ version, it would be a logical development to inherit the class of 
pseudo-random number generators TPSEMAR from a unique abstract class, and in this way 
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to define a universal interface for a library of tlie number generators. We intend to collect 
a library of a few random number generators with a universal interface (or find one) for 
the use in Foam and applications based on it. 

Concerning the version of Foam adapted to parallel processing, as in Refs. P,|8HT0|, 
we do not have plans in this direction in the immediate future. Here, we have of course 
in mind the use of the true CPU parallelism in the foam of cells build-up. One has to 
remember that in the high energy physics applications, which are our main objective, 
the foam of cell build-up will always be a tiny fraction of the total CPU time. The 
main fraction will be the subsequent MC simulation in which, as the vast experience with 
PC-farms in CERN and FNAL shows, one may organize the MC simulation with a low 
level of parallelism, with many simulators started with different random seeds running 
in parallel but not communicating - another specialized job combines all the results at 
the very end of the run. However, the first practical examples of true parallelism in 
the massive MC simulation for the purpose of the high energy physics experiments has 
appeared recently and is used in the LEP2 data analysis. 

The main algorithm of Foam is quite stable at this stage, and the main emphasis in 
its future development will be on making it more user friendly and better adapted to the 
use as a part of bigger MC projects. In particular its provisions for multibranching will 
become more sophisticated, as more feedback comes from real-life applications. 



5 Usage of Foam 

In the following we provide basic information on how to use the Foam program in the user 
application. 

5.1 Foam distribution directory of the C-\ — |- version 

The Foam package is distributed together with the demonstration main programs and 
some utilities in the form of about 20 files in a single UNIX directory F0AM-export-v2 . 05. 
Demonstration runs can be executed using standard make commands as follows: 

make Demo -run 
make DemoPers 
make Demo -map 
make DemoNR-run 

The essential fragments of the output from running command "make Demo-run" are shown 
in Appendix B. The compilation and linking procedure is encoded in the file Makefile, 
and it has to be checked by the user whether it conforms to the local operating sys- 
tem. In particular, if ROOT is used, then certain paths and environmental variables in 
Makefile have to be adjusted. The use of ROOT is decided by the presence of the vari- 
able ROOT_DEF in the file ROOT_DEF.h. Without ROOT, the user should execute command 
"make DemoNR-run". 
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The essential part of Foam, that is class TFOAM, TFCELL and a few auxiliary classes, 
is located in the files TFOAM . cxx and TFOAM . h. This is the "core" of the Foam source. 
The source code of the other utility classes TFHST, TFMAXWT, TPSEMAR and TFDISTR are in 
separate files. The main programs are in files Demo . cxx and DemoPers . cxx. They should 
serve as a useful template for the user's own application based on Foam. 

There is also one Fortran?? source code, circe2.fQ, which contains a testing den- 
sity distribution (2- dimensional beamstrahlung spectrum p5| of the electron linear col- 
lider 1^61) used in TFDISTR. The Maikef ile provides, therefore, also a useful example of 
linking C++ and F77 codes. 

There are also two output files output-Demo . linux and output -DemoNR. linux, which 
the reader may use to check whether he is able to reproduce these benchmark output 
results. 



5.2 Simple example of an application 

A very simple example of the use of Foam may look as follows: 

// *** Initialization *** 
double MCwt; 

TFDISTR *Densityl = new TFDISTR (FunType) ; // Create integrand function 
TPSEMAR *PseRan = new TPSEMARO ; // Create random numb, generator 

TFOAM *FoamX = new TFOAM ( "FoamX" ) ; // Create Simulator 
FoamX->SetkDim( 3); // Set dimension, h-rect . 

FoamX->Initialize(PseRan, Densityl ); // Initialize simulator 
// *** MC Generation *** 

TFHST *hst_Wt = new TFHST (0 . , 1 . 25, 25); // Create weight histogram 
double *MCvect =new double [3] ; // Monte Carlo event 

forClong loop=0; loop<1000000; loop++)-[ 

MCwt = FoamX->MCgenerate (double *MCvect) ; // Generate MC event 
hst_Wt->Fill(MCwt,1.0) ; // Fill weight histogram 

} 

// *** Finalization *** 
double IntNorm, Errel; 

FoamX->Finalize( IntNorm, Errel); // Print statistics, get normalization 
double MCresult, MCerror, AveWt, WtMax, Sigma; 

FoamX->GetIntegMC( MCresult, MCerror); // get MC integral 

double eps = 0.0005; 

FoamX->GetWtParams (eps , AveWt, WtMax, Sigma); // get MC wt parameters 
hst_Wt->Print ; // Print weight histogram 

The user normally provides his own density distribution function belonging to the 
class that has to inherit from the following abstract class: 
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We thank Thorsten Ohl for providing us with a prehminary version of this code K 
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class TFOAM_INTEGRAND{ // Abstract class of distributios for Foam 
public : 
TFOAM.INTEGRANDO { }; 
virtual "TFOAM.INTEGRANDO { }; 
virtual double Density(int ndim, double*) = 0; 

}; 

In the above example the distribution *Densityl belongs to the class TFDISTR, which is 
provided in the Foam distribution directory. 



Par am. 


Value 


Meaning 


llUllli 


n* 
u 


UlilicllbiOil Ol lilt: blliipilCcll bUUbpciCc 


1.^ 1 1 1 Wl 

is.LJLlLL 


n* 
u 


L-'llliclloiUii Ui Llifc: liyjJcllcCLclllgUlcii oUUoJJclCfc; 


nCplls 


1 nnn* 


A/Ta'VTrmTm ■nmnnPT' or r'f^llt; 


nSampl 


200* 


No. of MC events in the cell MC exploration 


nBin 


8* 


No. of bins in edge-histogram in cell exploration 


OptRej 


0* 


OptRej = 0, weighted; =1, w = 1 MC events 


OptDrive 


2* 


Maximum weight reduction 




1 


or variance reduction 


OptPeek 


0* 


Next cell for split with maximum i?^ (PeekMax) 




1 


or randomly with probability ~ R'^ (PeekRan) 


OptEdge 


0* 


Vertices are NOT included in the cell MC exploration 




1 


or vertices are included in the cell MC exploration 


OptOrd 


0* 


Root cell is a hyperrectangle in simplical subspace 




1 


or root cell is a simplex in simplical subspace 


OptMCell 


1* 


Economic memory algorithm in hyperrectangular subspace is ON 







or economic memory algorithm in hyperrectangular subspace is OFF 


EvPerBin 


25* 


Maximum number of effective w = 1 events/bin 







or counting of number eff. events/bin is inactive 


Chat 


1* 


=0,1,2 is the "chat level" in the standard output 


MaxWtRej 


1.1* 


Maximum weight used to get w = 1 MC events 


Table 6: Fourteen principal configuration parameters and switches of the Foam program. The 



default values are marked with the star superscript. 



5.3 Configuring Foam 

Foam has fourteen principal configuration parameters plus parameters inhibiting and/or 
predefining the division geometry in the cell split. 

5.3.1 Principal configuration parameters 

All of the principal parameters listed in Table |^ are set to meaningful default values, 
hence the beginner may stay ignorant of their role for some time, and learn gradually how 
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to exploit them in order to improve the efficiency of Foam in the actual application. All 
these parameters are data members of the TFOAM class, see Table. If the user wants to 
redefine all of them, then the relevant piece of code will look as follows: 

FoamX->SetnDim( nDim) ; 
FoamX->SetkDim( kDim) ; 
FoamX->SetnCells( nCells) ; 
FoamX->SetnSampl ( nSampl) ; 
FoamX->SetnBin( nBin) ; 
FoamX->SetOptRe j ( OptRe j ) ; 
FoamX->SetOptDrive ( OptDrive) ; 
FoamX->SetOptPeek( OptPeek) ; 
FoamX->SetOptEdge( OptEdge) ; 
FoamX->SetOptOrd( OptOrd) ; 
FoamX->SetOptMCell( OptMCell) 
FoamX->SetEvPerBin( EvPerBin) 
FoamX->SetMaxWtRe j ( MaxWtRe j ) 
FoamX->SetChat ( Chat); 

In practical applications one will redefine only some of them. The minimum requirement 
is that the user sets a non-zero value of nDim or kDim such that the total dimension 
nDim+kDim is a non-zero positive integer. 



5.3.2 Inhibiting cell division in certain directions 

If the user of Foam decides to inhibit the division in certain variable in the hyper rectangular 
subspace, this can be done with the method SetInhiDiv(int iDim, int InhiDiv) of 
the class TFOAM, where iDim is the dimension index for which inhibition is done and 
InhiDiv is the inhibition tag. This method should be used before invoking Initialize, 
after setting nDim and/or kDim. The relevant code may look as follows: 

FoamX->SetInhiDiv(0 , 1); //Inhibit division of x_l 
FoamX->SetInhiDiv(l , 1); //Inhibit division of x_2 

The allowed values are InhiDiv=0,l and the default value is InhiDiv=0. Note that the 
numbering of dimensions with iDim starts from zero and variables of the hyperrectangular 
subspace always come ffist, before the simplical ones. 



5.3.3 Setting predefined cell division geometry 

The user may predefine divisions of the root cell in certain variables in the hyperrectan- 
gular subspace using the method SetXdivPRD(int iDim, int len, double xDiv[]). 
The relevant piece of the user code may look as follows: 

double xDiv[3] ; 

xDiv[0]=0.30; xDiv [1] =0 . 40 ; xDiv [2] =0 . 65 ; 
FoaiiiX->SetXdivPRD(0, 3, xDiv) ; 
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Again, this should be done before invoking Initialize, after setting nDim and/or kDim. 



5.4 Persistency 



Persistency of the Foam classes is arranged using "default streamers" of the ROOT pG 
package. Writing a TFOAM class object into the disk file rmain.root can be done with the 
single Write as follows: 

TFile RootFileC "rmain.root" , "RECREATE" , "histograms") ; 

FoamX->Write("FoamX") ; //Writing Foam on the disk 

RootFile.WriteO ; 
RootFile . Close ; 

The instruction FoamX->Write("FoamX") can be put at any place of the code after the 



instruction FoamX->Initialize( . . . ), see example of the user code shown in Section 

Next, in another program, the TFOAM class object can be read from the disk file 
rmain.root as follows: 

TFile f ileA("rmain.root") ; // connect disk file 
f ileA.cdO ; 

fileA.lsO; // optional printout 

fileA.MapO; // optional printout 

f ileA . ShowStreamerlnf ; // optional printout 

f ileA.GetListOfKeys()->Print() ; // optional printout 
TFOAM *FoamX = (TFOAM*)f ileA.Get("FoamX") ; // find object 
FoamX->LinkCells ; // restore pointers of the binary tree of cells 
FoamX->CheckAll(l) ; // optional x-check of pointers 

and at this point the FoamX object is ready to generate MC events, as in the MC generation 
part of the code shown in Section [5l^ . 



5.5 Fortran77 versions 

The distribution directory FoamF77-2 . 02-export contains the README file, two demon- 
stration main programs DemoFoam.f and DemoMCell.f to be compiled and run with the 
help of commands 

make DemoFoam 
make DemoMCell 

encoded in the Makefile. The outputs from the above runs can be compared with the 
benchmark outputs output-DemoFoam-linux and output-DemoMCell-linux. 

The basic Foam source files are: FoamA.f with header file FoamA.h and MCellA.f with 
header file MCellA.h. 
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For the description of the input (configuration) parameters, see comments in FoamA . f 
and MCellA.f respectively. The names of the configuration variables are the same as in 
the C++ version, except nCells, which is here renamed to nBuff . Their values and the 
meaning are the same. 

Demonstration main programs DemoFoam.f and DemoMCell.f can serve as templates 
for the user application programs. 

The testing main program uses the histogramming package GLK of the KKMC pro- 
gram [|l^], which the user may replace with any other histogramming package. 




Figure 15: Testing distribution Pcamei{x) of Ref. Q in two dimensions. 



6 Numerical studies and example applications 

In the following subsection we examine MC efficiency of the Foam in a series of numerical 
exercises. In some of them we shall also show examples of the Foam application with the 
distributions relevant to everyday practice in high energy physics. 

6.1 Dependence of the Foam efficiency on the configuration pa- 
rameters 

As a first numerical exercise, we examine the dependence of the Foam efficiency on the 
most important (input) configuration parameters, including the dimension of the space. 

In Table we collect results from many MC runs for various dimensions, numbers of 
cells and numbers of MC events in single cell exploration, varying also the type of the 
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nDim 


kDim 


nCalls 


nCells 


nSampl 




a/{w) 


^statist. R 


R 





1 


203719 


1000 


333 


0.99148 


0.014585 


1.031e-05 


0.99999782 





1 


206192 


1000 


1000 


0.99147 


0.014752 


1.043e-05 


0.99999962 





1 


206192 


1000 


3333 


0.99147 


0.014752 


1.043e-05 


0.99999962 





1 


206192 


1000 


10000 


0.99147 


0.014752 


1.043e-05 


0.99999962 





1 


206192 


1000 


33333 


0.99147 


0.014752 


1.043e-05 


0.99999962 





3 


275421 


1000 


333 


0.50538 


0.54088 


0.000382 


0.99986832 





3 


435112 


1000 


1000 


0.50886 


0.54033 


0.000382 


1.00017104 





3 


663493 


1000 


3333 


0.49922 


0.55721 


0.000394 


0.99934710 





3 


834094 


1000 


10000 


0.50359 


0.54674 


0.000386 


1.00056316 





3 


1015157 


1000 


33333 


0.51091 


0.54035 


0.000382 


0.99983999 





3 


2312759 


10000 


333 


0.72346 


0.27758 


0.000196 


0.99977045 





3 


2675820 


10000 


1000 


0.72677 


0.27504 


0.000194 


0.99995080 





3 


3054404 


10000 


3333 


0.72199 


0.27710 


0.000195 


1.00013270 





3 


3333479 


10000 


10000 


0.72200 


0.27720 


0.000196 


1.00008994 





3 


3575366 


10000 


33333 


0.72243 


0.27786 


0.000196 


0.99997875 





4 


3825046 


10000 


1000 


0.50363 


0.51168 


0.000361 


1.00013082 





4 


6559430 


10000 


10000 


0.50297 


0.51001 


0.000360 


0.99960319 


2 


2 


4493961 


10000 


1000 


0.43076 


0.63185 


0.000446 


1.00072564 


2 


2 


9374351 


10000 


10000 


0.44922 


0.60669 


0.000429 


1.00013171 


4 





6642202 


10000 


1000 


0.21029 


1.19420 


0.000844 


1.00072248 


4 





12337748 


10000 


10000 


0.20817 


1.20067 


0.000849 


1.00020405 





6 


2311881 


1000 


3333 


0.04199 


2.12091 


0.001499 


0.99856206 





6 


5542146 


1000 


10000 


0.03847 


2.38588 


0.001687 


0.99912901 





6 


12844256 


1000 


33333 


0.03279 


2.61028 


0.001845 


0.99799089 





6 


12737314 


10000 


3333 


0.15385 


1.15211 


0.000814 


1.00039754 





6 


24134694 


10000 


10000 


0.15313 


1.19596 


0.000845 


0.99945766 





6 


42827237 


10000 


33333 


0.14168 


1.22627 


0.000867 


0.99954178 





6 


42808972 


100000 


1000 


0.30910 


0.71250 


0.000503 


0.99972833 





6 


61803017 


100000 


3333 


0.30805 


0.71462 


0.000505 


1.00002674 





6 


92531875 


100000 


10000 


0.30905 


0.71423 


0.000505 


0.99985093 





9 


78325890 


100000 


1000 


0.03718 


1.64608 


0.001163 


0.99367339 





9 


167710365 


100000 


3333 


0.05247 


1.73063 


0.001223 


1.00109792 





9 


353943409 


100000 


10000 


0.05196 


1.80538 


0.001276 


1.00196909 





9 


272162624 


400000 


1000 


0.08490 


1.30193 


0.000920 


1.00065580 





9 


495260998 


400000 


3333 


0.09174 


1.35307 


0.000956 


0.99884358 





9 


924011087 


400000 


10000 


0.08853 


1.38579 


0.000979 


1.00052122 





12 


261911066 


100000 


3333 





5.83954 


0.004129 


0.97304842 





12 


671460574 


100000 


10000 


0.00640 


3.85823 


0.002728 


0.98878698 





12 


913072065 


400000 


3333 


0.01285 


2.73991 


0.001937 


0.98688299 





12 


2117963809 


400000 


10000 


0.01235 


2.92642 


0.002069 


0.99301117 



Table 7: Numerical results of Foam with the maximum weight reduction. Variable nCalls is 
the total number of the function calls in the foam build-up. 
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nDim 


kDim 


n Calls 


nCells 


nSampl 




a/{w) 


^statist. R 


R 





4 


3855289 


10000 


1000 


0.27659 


0.31944 


0.000225 


1.00025027 





4 


7760907 


10000 


10000 


0.30313 


0.31483 


0.000222 


0.99978048 


2 


2 


4589024 


10000 


1000 


0.23086 


0.38050 


0.000269 


0.99967167 


2 


2 


8696153 


10000 


10000 


0.24696 


0.37153 


0.000262 


0.99959278 


4 





6157799 


10000 


1000 


0.08498 


0.92314 


0.000652 


1.00006553 


4 





10547749 


10000 


10000 


0.09881 


0.89859 


0.000635 


1.00024727 



Table 8: Numerical results of Foam with the variance reduction. 



cells. We always use the same test distribution pcamei{x) as in Ref. 0, which features 
two relatively narrow gaussian peaks placed on the diagonal. The 2-dimensional version 
of this distribution is shown in Fig. We have used the non-default values nBin=4 and 
EvPerBin=50 and the default values for the other configuration parameters. In this table 
all tests were done for the default option of reduction of the maximum weight, 0ptDrive=2 
(for results with OptDrive=l, see next table). The efficiency w^ax/ ('^) is calculated using 
maximum weight w^^^ defined as inQ Ref. [IH for £ = 0.0005. The maximum weight w^ax 
is calculated with the help of the small auxiliary class TFMAXWT. 

In Table |^ the efficiency of the MC run measured in terms of Wmax/ {w) and a/ {w). 
The value of the integral R ± ^statist.Ri shown in the last four columns, was obtained 
from the MC run in which the total number of MC events was 2 x 10^. The value of the 
integral R is well known, it is equal to 1, within 10~^. 

The following observations based, on the results of Table |^, can be made: 

• Looking at the results, for total dimension n = 4 we see that the hyperrectangular 
cells clearly provide better MC efficiency than simplical ones. All other results are 
for hyperrectangular cells. 

• All of the results support the observation that the MC efficiency depends critically 
on the number of cells. In particular, see results for n = 6, the increase of nSamp 
(No. of MC events in cell exploration) beyond a certain value does not help at all. 

• In the case of the very inefficient Foam, see n = 12 with cr/{w) ~ 6, the estimate of 
the MC statistical error can be misleading. We see an indication that one should 
not trust runs with a/ {w) > 3. 

• For this particular testing function the dimension n = 12 requires a minimum of 
400k cells and resulting MC efficiency of order of 1% is barely acceptable^. 

In Table § we repeat the exercise of Table |^ for the option of the variance reduction 
OptDrive=l at four dimensions. As compared to Table ^ we see a net improvement in 
the variance and deterioration of the w'^^^. This agrees with the expectations. 

^^The e-dependent maximum weight is defined such that events with w > w^ax contribute an £-fraction 
to the total integral. It is numerically more stable in the numerical evaluation than the one defined as 
the largest weight in the MC run. 

^^However, we still get the correct value of the integral within 0.2%. 
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Functions at 2-dimens. 


Foam 1.01 


Simpl. 


H-Rect. 


VEGAS 


Pa{x) (diagonal ridge) 


0.93 


0.93 


0.86 


0.03 


Pb{x) (circular ridge) 


0.82 


0.82 


0.82 


0.16 


Pc{x) (edge of square) 


0.57 


1.00 


1.00 


0.53 


Functions at 3-dimens. 


Foam 1.01 


Simpl. 


H-Rect. 


VEGAS 


Pa{x) (thin diagonal) 


0.67 


0.74 


0.66 


0.002 


Pb{x) (thin sphere) 


0.36 


0.47 


0.53 


0.11 


Pc{x) (surface of cube) 


0.37 


0.95 


1.00 


0.30 



Table 9: Efficiencies (w)/Wmax ^ = 0.0005. Functions pxix) are the same as in Ref. [p. 
Results from Foam are for 5000 cells (2500 active cells) and cell exploration is done for a 
modest 200 MC events/cell. 





k 


n 


Nc 


Ns 


Nb 


bin 


nCall 


M 


(w) 

"'max 


R±AR 


AR/R 


1 


2 





IK 


IK 


4 


25 


900K 


1168.8 


0.0 


5.40726±1.99871 


0.36963 


1 

2 






2 
2 


IK 
IK 


IK 
IK 


4 
4 


25 
25 


169K 
215K 


0.2272 
0.2962 


0.6149 
0.7754 


3.14121±0.00022 
3.14118±0.00029 


7.1-10-^ 
9.3-10-5 


1 
1 






2 
2 


5K 
lOK 


IK 
IK 


4 
4 


25 
25 


656K 
1174K 


0.0639 
0.0487 


0.8421 
0.8877 


3.14159±0.00006 
3.14156±0.00005 


2.0-10-5 
1.5-10-5 


1 
1 






2 
2 


IK 

5K 


lOK 
lOK 


4 
4 


25 
25 


849K 
1457K 


0.1479 
0.0606 


0.5920 
0.8354 


3.14118±0.00014 
3.14150±0.00006 


4.6-10-5 
1.9-10-5 


1 
1 






2 
2 


IK 
IK 


2K 
8K 


8 
8 


25 
100 


621K 
1671K 


0.0606 
0.1048 


0.8354 
0.6652 


3.14195±0.00026 
3.14168±0.00010 


8.4-10-5 
3.3-10-5 



Table 10: Numerical results of Foam for 2-dimensional distribution of Eq. (|^) for p = 10 ^. 
Variation of the configuration parameters: k =kDim, n =nDim, =nCells (No. of function 
calls), Nb =nBirL, ^^^^ =EvPerBin. In the first column we mark the type of the weight 
optimization OptDrive=l ,2, for variance or maximum weight reduction. The value of the 
integral R and its statistical error AR are from MC run of Nmc = 10"^ events, w^ax 
e = 0.0005. nCalls is the total number of the function calls in the foam build-up. 

6.2 Comparison with Foam 1.x and classic VEGAS 

In Table |^ we update the comparison of Foam and VEGAS of Ref. [|l| , adding results for the 
new hyp err ect angular option. The simplical results are now clearly improved with respect 
to Ref. |1|], because of the better cell division algorithm. Generally, the hyperrectangular 
cell mode provides as good an efficiency as the simplical one. However, one should keep in 
mind that Foam with hyperrectangular cells is a factor of 2 or more faster in the execution. 
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6.3 Example of sharply peaked distribution 



In Table ^ we examine the dependence of the MC efficiency/error on the various input 
configuration parameters of Foam. All these numerical results are for the distribution 

which, for /i = 10~^, has a very sharp ridge across the diagonal Xi + X2 = 1. This 
distribution is taken from Refs. |]I2|,|I^ and is related to the photon distribution at high 
energy electron-positron colliders. 




Figure 16: Rectangular and triangular foam of 1000 cells for the distribution of Eq. (p^). 

What can we learn from the results in Table |10]? First of all, in the first line, we see 
a spectacular failure of Foam with rectangular cells0. The value of the integral is wrong 
by a factor of 2 and the statistical error is underestimated. This illustrates the problem 
of the the lack of "angular mobility" of the rectangular cells discussed in Section |2.3| ; 
rectangles are unable to align with the singularity along the diagonal. 

This we illustrate in the left plot of Fig. |T6|, for rectangular 1000 cells, where we clearly 
see "blind spots" . In the right-hand plot the triangular cells in the foam are clearly aligned 



with the diagonal ridge. In the rows 2 and 3 of Table |10] we see the reasonable numerical 
results for the triangular foam, which are for the maximum weight reduction and variance 
reduction options respectively; the other configuration parameters are rather close to the 
default ones. In rows 4 and 5 we are playing with the increase of the cell number and in 
the rows 6 and 7 with the number of MC events used in the cell exploration. Finally in 
rows 8 and 9 we change the binning of the histograms used in the MC cell exploration. 
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Setting the number of the MC events in a single rectangle exploration to 10^ cures the problem partly. 
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As we see, the most profitable in terms of MC efficiency /precision is the increase of the 
number of cells however, adjusting other parameters can also help. In all cases we display 
variable nCalls, the number of calls of the density distribution function during the foam 
build-up. In the best result of line 5 with 10000 triangular cells we have obtained a 5-digit 
precision for about 10^ function calls^. 

Summarizing, we see from the above exercise that the user of Foam has a possibility to 
adjust several configuration parameters, so that the MC efficiency for a given distribution 
is improved quite significantly. 



MAPPING 


k 


n 


Nc 


Ns 


Nb 


We// 

bin 


nCall 




(w) 

•^max 


AR/R 


SIZE 


(al) OFF 
(a2) ON 


8 
9 






1 
1 


IK 
IK 


4 
4 


25 
25 


518 
218 


2.1555 
1.1391 


0.038 
0.115 


0.00481 
0.00254 


15KB 
15KB 


(bl) OFF 
(b2) ON 


8 
9 






20 
20 


IK 
IK 


4 

4 


25 
25 


5767 
3548 


1.1847 
0.7626 


0.130 
0.206 


0.00264 
0.00170 


15KB 
15KB 


(cO) OFF 
(cl) OFF 
(c2) ON 



8 
9 


8 




1000 
1000 
1000 


IK 
IK 
IK 


4 
4 
4 


25 
25 
25 


487K 
145K 
125K 


1.6603 
0.5298 
0.7626 


0.085 
0.359 
0.394 


0.00371 
0.00118 
0.00104 


54KB 
53KB 
53KB 


(dl) OFF 
(d2) ON 


8 
9 






5000 
5000 


IK 
IK 


4 
4 


25 
25 


528K 
596K 


0.4330 
0.4037 


0.438 
0.467 


0.00096 
0.00090 


209KB 
209KB 



Table 11: Numerical results from Foam simulation/integration for the decay process r — > 
z/TT^TT+TT^, according to the matrix element squared of the TAUOLA MC event generator P5| , P^ |. 
All MC averages are for 200K events generated after Foam initialization with the configuration 
parameters given in the table. The notation and the meaning of the quantities are the same 
as in previous tables. The size of the ROOT disk file, in which the Foam object was written is 
indicated in the last column. 



6.4 Decay of the r lepton into 3 pions 

In Table |Tl| we collect numerical results for an example of the application of Foam to the 
very practical problem of the MC simulation of the decay r —>■ wk'tt^tt" , according to 
the distribution (matrix element squared) taken from the MC program TAUOLA p8| , |29[] . 

The amplitude of the decay process contains two distinct parts due to two Feynman 
diagrams, see Fig. |1^, which have peaks due to the ai resonance and the p resonance. 
There are two peaks due to p resonances that partly overlap in the Lorentz invariant phase- 
space, such that the actual shape of the differential distribution is rather complicated. 
We took for this exercise subroutine DPHTRE of TAUOLA, in which nine random numbers^ 
are replaced by the nine variables of Foam. The 4-particle phase-space is 8-dimensional. 
The ninth variable is due to two branches in the phase-space parametrization of TAUOLA, 
and in case of Foam as well (the method is similar to that of Section |2.11| ). In cases of "no 

^°In Ref. the same precision for the same function was attained for about 10^ function calls. 
^^Including two random numbers of the subroutine SPHERA and two (Euler) angles corresponding to 
the overall rotation of the entire event. 
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Figure 17: Feynman diagrams for r decay into 3 pions. 




Figure 18: The distribution of r — > uStt decay as a function of xi and X2 (without mapping), fixing 
the other variables to some values. 

mapping and no multibranching" we are back to eight dimensions. The variables xi and 
X2 of Foam represent (up to a hnear transformation) the two effective masses of the Sir 
and 27C systems. The next four variables Xj, i = 3, 4, 5, 6 are the polar variables cos^ and 
of the pions in the rest frame of the Svr and 2tt systems - this is a completely standard 



52 



phase-space parametrization, see Ref. [|^, and also Ref. The variables xj and xs are 



reserved for the overall rotation, and the last one, xg, is mapped into branch index j = 1,2, 



see Section Variables Xj and Xg are inhibited (no cell-split in them), because the 

distribution does not depend on them. Variable Xg (if present) has a predefined division 
value equal to 0.5 and is inhibited for the division (see Section In Fig. |18| we show 

the decay distribution as a function of xi and X2 (without mapping), fixing the other 
variables to some values. The distribution is clearly a non-trivial one. 

In Table we show the MC efficiency of Foam for a gradually increasing number of 
hyperrect angular cells, with the mapping compensating for the Breit-Wigner peaks of the 
ai and p resonances (as in Ref. ||29|) and without. One example with simplical cells is 
also included. 

The most striking result in Table |TT| is the comparison of lines (a2) and (bl): the Foam 
algorithm with only 20 cells is performing equally well as the doubly-branched mapping 
compensating for the resonance peaks of the ai and p. When going to higher number of 
cells, the MC efficiency in the cases with and without mapping becomes almost the same. 
This is expected, since Foam also does the mapping compensating for the resonances on its 
own. From row (cO) we see also that the simplical mode of Foam is clearly underperforming. 
We think that Foam configured with 1000 cells, see row (cl) in Table |Tl|, is an economic 
solution to this particular MC problem^. (Mapping due to resonances is not really 
necessary in this case.) 



^^The net profit with respect to TAUOLA ^ would be a three times faster program, and more 
importantly, a significantly simpler code. 
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Figure 19: Foam of rectangular cells for the electron-positron beamstrahlung spectrum of 
circe2 [||] with 500 cells. 



6.5 Beamstrahlung spectrum 

In Fig. we show the foam of rectangular cells for the 2-dimensional beamstrahlung 
spectrum D{zi, Z2) of the electron-positron collider at 500 GeV, encoded in the pro- 
gram circe2 of Refs. |2^, It should be stressed that this spectrum is not known 



analytically but rather through a numerical fit to results of the machine simulation or 
(in the future) from an experiment. In order to avoid (integrable) infinite singularities at 
Zi = I'm D{zi, Z2) we use the variables tj = (1 — i = 1, 2 in Fig. 19. For this exercise 
we used foam of 500 cells getting the MC efficiency cr j {w) = 0.41 and {w)/w^^^ = 0.64 
(for e = 0.0005); enough for practical application (can be easily improved by adding more 
cells) . 

Generating D{zi, Z2) is not really such a very much important and difficult problem. 
A more interesting problem is to generate the distribution D{zi, Z2)(t{sziZ2), where a{s) 
is the cross section of some physics process, which may have a strong singularity of its 
own, such as a resonance or threshold factor. Such a problem was already treated with 



the help of the Foam program in KKMC |jT8[ event generator and the study of Ref. |^| 
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7 Conclusions 

The author hopes that this new adaptive tool for constructing efficient MC programs will 
find its way to many applications in high energy physics and beyond. The main points 
on the new Foam algorithm and the program are the following: 

• Foam is a versatile adaptive, general-purpose Monte Carlo simulator. 

• The Foam algorithm is based on the cellular division of the integration domain. 

• The geometry of the "foam of cells" is rather simple, cells of simplical and/or hy- 
perrectangular shape are constructed in the process of a binary split. 

• It works in principle for arbitrary distribution - no assumption of factorizability as 
in VEGAS of Ref. §. 

• Foam is reducing maximum weight of the weight distribution, it can therefore provide 
unweighted events. The variance reduction, useful for the integration and generating 
weighted events, is also available. 

• Encoding cells with memory- efficient methods allows up to ~ 10^ cells to be built 
in the computer memory of a typical desktop computer. 

• The rules for choosing the next cell for division and the division geometry are based 
on rather sophisticated analysis of the projection on the cell edges. This costs CPU 
time, which becomes the main barrier towards higher MC efficiency. 

• Foam can deal efficiently with rather strongly peaked distributions, up to relatively 
high dimensions, ~ 12 dimensions, with today's desktop computers. 
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A Variance optimization 



Suppose we have already constructed the cells c<Ji, c<Ji, . . . , cj^v and within each cell we 
have defined the function p'{x) constant over the cell, p'{x) = p'j. The integral R'j = 
p'{x)dx"' = p'jVi is known, because the volume of the cell Vj is known. The function 
p{x) is not, in general, constant over the cell and the weight w = p{x)/p'{x) is used to 
determine its integral in the usual way: R = R'{w)p' = ^^jiw) p'^Ri, where the average 



(^)p'i = 4r / p'i{x)a{x)dx'' 



is defined for the J-th cell alone. 

The question is now the following: preserving the geometry of the cells, can we get 
smaller variance by simply changing the probabilities of the generation of the cells? Rescal- 
ing p'j —>■ p'l = Xip'j does not affect the integral R, because the change of the normalization 
of R' and {'w)p/ is cancelling. It is convenient to assume that the above rescaling preserves 
R' = R" and the total average weight {w)pi = {w)p", that is A/ obey the constraint 
^jR'i = = const. With the above constraint in mind, we now ask: for what 

values of A/ is the dispersion of the total weight cr^ = (w^)p" — minimal?]^ Since 

by construction {w)p>r is independent of A/, we may look only for a minimum of {w'^)p". 
With the standard methods we get a (local) minimum condition: 



dXiV ) dX 



where A is the Lagrange multiplier. The solution of the minimum condition 



(27) 



§(^^^^-Ai?; = (28) 



is simply A/ ~ const x \/ {w^)p' , or more precisely 



A/ = _ (29) 



The value of {w^)ram = (^^)p" calculated at the minimum, that is for A/ of Eq. (ppj), is 
also rather simple 



This problem was, of course, often considered in the past, see for instance Refs. |,20|. We outhne 
here the sohition for the sake of completeness and convenience of the reader. 
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Let us note that the functional a/ (tf^)min, which we intend to minimize in the process 
of the evolution of the foam (cell split) is a simple sum of contributions from all cells. 
Consequently, when working out details of the split of a given cell we may calculate the 
gain in terms of the total dispersion independently of other cells, see also below. This 
very convenient feature is exploited in the algorithm of the foam build-up. 

Adopting (temporarily) the following normalization conventions: p'j = 1, w = p{x) 
and R'l = Vj, Xj ~ we get 



P"i = VWi, Ri = Vi^/Wi. (31) 

where the average (...)/ is understood as defined for points uniformly distributed within 
the J-th cell. 

In the Foam algorithm we do not go, of course, through a judicious adjustment of the 
relative importance of the cells, which minimizes the variance as described above; but 
instead, we simply declare that the distribution p\x) is defined by the optimum solution 
of Eq. (prp. Once we have done that, for the MC weight w defined with respect to such a 



new p\x), we find out that {w )i = 1 for each cell and, also for all cells, {w ) = 1. What 
is then minimized in the process of the cell division is the ratio of the variance to the 
average weightQ 



This quantity is not so convenient to optimize in the process of the cell split (foam 
evolution) and we rather chose to minimize a closely related 'linearized" quantity 



■loss 




L j = ^/(/R)7 - {w)i) = j pioss{x)dx'^ = R'-R, 



(33) 

which is a sum of contributions from all cells and is a monotonous ascending function of 
cr/ {w). In the process of the cell division uj ujia®^ib, we decrease a / {w) step by step, 
by playing with the geometry of the cell split, so that the gain in the total Rioss = R' — R 
due to a given cell split is as big as possible. It is convenient that the contributions from 
the other cells to Rioss are unchanged. In this way every cell split will lead to a smaller 
and smaller a/ (w) in the final MC run. 



We exploit here the relation (w) — R/R', and we should keep in mind that R is constant during the 
variance minimization. 
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B Output of the demonstration program in CH — \- 



T7 
r 

p 








r 








TFOAM: : Initialize ****** 


r 








p 










FoamX 


p 


V CX o X Ull 






2. 


05 = Release date: 2002.02.13 


p 













TT 
r 






1 




Dimension of the hyper-cubical sub-space 


p 






1000 




Requested number of Cells (half of them active) 


F 


nSampl 


= 


500 




No of MC events in exploration of a cell cell 


F 


nBin 




8 




No of bins in histograms, MC exploration of cell 


F 


EvPerBin 




25 




Maximum No ef f ective_events/bin, MC exploration 


F 


OptDrive 




2 




Type of Driver =1,2 for Sigma, WtMax 


F 


OptEdge 









Decides whether vertices are included in the MC 


F 


OptPeek 









Type of the cell Peek =0,1 for maximum, random 


F 


OptOrd 









Root cell hyp-cub. or simplex, =0,1 


F 


OptMCell 




1 




MegaCell option, slim memory for hyp-cubes 


F 


OptDebug 




1 




Additional debug histogram, SetDirectory (1) 


F 


OptCulst 




1 




Numbering of dimensions starts with h-cubic 


F 


OptRej 









MC rejection on/off for OptRej =0,1 


F 
F 


MaxWtRej 








2 = Maximum wt in rejection for wt=l evts 



FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF 



F 

F 
F 
F 
F 
F 
F 



nCalls = 
XPrime = 

XDiver = 
MCresult = 



*** TFOAM: : Initialize FINISHED! ! ! *** 



246799 = 
1 . 2748942 
0.27441356 
1.0004806 = 



Total number of function calls 
Primary total integral 
Driver total integral 
Estimate of the true MC Integral 



F 
F 
F 
F 
F 
F 
F 
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FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF 

F F 

F ****** TFOAM: : Finalize ****** p 

F NevGen = 2000 = Number of generated events in the MC generation F 

F LastVe = -1 = Number of vertices (only for simplical option) F 

F nCalls = 248799 = Total number of function calls F 

F F 

F AveWt = 0.78836949 = Average MC weight F 

F WtMin = 1 . 3022085e-06 = Minimum MC weight (absolute) F 

F WtMax = 1.0227719 = Maximum MC weight (absolute) F 

F F 

F XPrime = 1.2748942 = Primary total integral, R_prime F 

F XDiver = 0.27441356 = Driver total integral, R_loss F 

F F 

F IntMC = 1.0050877 +- 0.0052357675 = Result of the MC Integral F 

F MCerelat = 0.0052092644 = Relative error of the MC intgral F 

F <w>/WtMax = 0.80037512 = MC efficiency, acceptance rate F 

F Sigma/<w> = 0.23296539 = MC efficiency, variance/ave_wt F 

F WtMax = 0.985 = WtMax (esp= 0.0005) F 

F Sigma = 0.1836628 = variance of MC weight F 

F F 
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